home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / filter.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-19  |  12.0 KB  |  583 lines
  1. /*
  2.  *        Filtered Image Rescaling
  3.  *
  4.  *          by Dale Schumacher
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include <string.h>
  9. #include <malloc.h>
  10. #include <math.h>
  11. #include "GraphicsGems.h"
  12.  
  13. static char    _Program[] = "fzoom";
  14. static char    _Version[] = "0.20";
  15. static char    _Copyright[] = "Public Domain 1991 by Dale Schumacher";
  16.  
  17. #ifndef EXIT_SUCCESS
  18. #define    EXIT_SUCCESS    (0)
  19. #define    EXIT_FAILURE    (1)
  20. #endif
  21.  
  22. typedef    unsigned char    Pixel;
  23.  
  24. typedef struct {
  25.     int    xsize;        /* horizontal size of the image in Pixels */
  26.     int    ysize;        /* vertical size of the image in Pixels */
  27.     Pixel *    data;        /* pointer to first scanline of image */
  28.     int    span;        /* byte offset between two scanlines */
  29. } Image;
  30.  
  31. #define    WHITE_PIXEL    (255)
  32. #define    BLACK_PIXEL    (0)
  33.  
  34. /*
  35.  *    generic image access and i/o support routines
  36.  */
  37.  
  38. static char *
  39. next_token(f)
  40. FILE *f;
  41. {
  42.     static char delim[] = " \t\r\n";
  43.     static char *t = NULL;
  44.     static char lnbuf[256];
  45.     char *p;
  46.  
  47.     while(t == NULL) {            /* nothing in the buffer */
  48.         if(fgets(lnbuf, sizeof(lnbuf), f)) {    /* read a line */
  49.             if(p = strchr(lnbuf, '#')) {    /* clip any comment */
  50.                 *p = '\0';
  51.             }
  52.             t = strtok(lnbuf, delim);    /* get first token */
  53.         } else {
  54.             return(NULL);
  55.         }
  56.     }
  57.     p = t;
  58.     t = strtok((char *)NULL, delim);        /* get next token */
  59.     return(p);
  60. }
  61.  
  62. Image *
  63. new_image(xsize, ysize)    /* create a blank image */
  64. int xsize, ysize;
  65. {
  66.     Image *image;
  67.  
  68.     if((image = (Image *)malloc(sizeof(Image)))
  69.     && (image->data = (Pixel *)calloc(ysize, xsize))) {
  70.         image->xsize = xsize;
  71.         image->ysize = ysize;
  72.         image->span = xsize;
  73.     }
  74.     return(image);
  75. }
  76.  
  77. void
  78. free_image(image)
  79. Image *image;
  80. {
  81.     free(image->data);
  82.     free(image);
  83. }
  84.  
  85. Image *
  86. load_image(f)        /* read image from file */
  87. FILE *f;
  88. {
  89.     char *p;
  90.     int width, height;
  91.     Image *image;
  92.  
  93.     if(((p = next_token(f)) && (strcmp(p, "Bm") == 0))
  94.     && ((p = next_token(f)) && ((width = atoi(p)) > 0))
  95.     && ((p = next_token(f)) && ((height = atoi(p)) > 0))
  96.     && ((p = next_token(f)) && (strcmp(p, "8") == 0))
  97.     && (image = new_image(width, height))
  98.     && (fread(image->data, width, height, f) == height)) {
  99.         return(image);        /* load successful */
  100.     } else {
  101.         return(NULL);        /* load failed */
  102.     }
  103. }
  104.  
  105. int
  106. save_image(f, image)    /* write image to file */
  107. FILE *f;
  108. Image *image;
  109. {
  110.     fprintf(f, "Bm # PXM 8-bit greyscale image\n");
  111.     fprintf(f, "%d %d 8 # width height depth\n",
  112.         image->xsize, image->ysize);
  113.     if(fwrite(image->data, image->xsize, image->ysize, f) == image->ysize) {
  114.         return(0);        /* save successful */
  115.     } else {
  116.         return(-1);        /* save failed */
  117.     }
  118. }
  119.  
  120. Pixel
  121. get_pixel(image, x, y)
  122. Image *image;
  123. int x, y;
  124. {
  125.     static Image *im = NULL;
  126.     static int yy = -1;
  127.     static Pixel *p = NULL;
  128.  
  129.     if((x < 0) || (x >= image->xsize) || (y < 0) || (y >= image->ysize)) {
  130.         return(0);
  131.     }
  132.     if((im != image) || (yy != y)) {
  133.         im = image;
  134.         yy = y;
  135.         p = image->data + (y * image->span);
  136.     }
  137.     return(p[x]);
  138. }
  139.  
  140. void
  141. get_row(row, image, y)
  142. Pixel *row;
  143. Image *image;
  144. int y;
  145. {
  146.     if((y < 0) || (y >= image->ysize)) {
  147.         return;
  148.     }
  149.     memcpy(row,
  150.         image->data + (y * image->span),
  151.         (sizeof(Pixel) * image->xsize));
  152. }
  153.  
  154. void
  155. get_column(column, image, x)
  156. Pixel *column;
  157. Image *image;
  158. int x;
  159. {
  160.     int i, d;
  161.     Pixel *p;
  162.  
  163.     if((x < 0) || (x >= image->xsize)) {
  164.         return;
  165.     }
  166.     d = image->span;
  167.     for(i = image->ysize, p = image->data + x; i-- > 0; p += d) {
  168.         *column++ = *p;
  169.     }
  170. }
  171.  
  172. Pixel
  173. put_pixel(image, x, y, data)
  174. Image *image;
  175. int x, y;
  176. Pixel data;
  177. {
  178.     static Image *im = NULL;
  179.     static int yy = -1;
  180.     static Pixel *p = NULL;
  181.  
  182.     if((x < 0) || (x >= image->xsize) || (y < 0) || (y >= image->ysize)) {
  183.         return(0);
  184.     }
  185.     if((im != image) || (yy != y)) {
  186.         im = image;
  187.         yy = y;
  188.         p = image->data + (y * image->span);
  189.     }
  190.     return(p[x] = data);
  191. }
  192.  
  193.  
  194. /*
  195.  *    filter function definitions
  196.  */
  197.  
  198. #define    filter_support        (1.0)
  199.  
  200. double
  201. filter(t)
  202. double t;
  203. {
  204.     /* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */
  205.     if(t < 0.0) t = -t;
  206.     if(t < 1.0) return((2.0 * t - 3.0) * t * t + 1.0);
  207.     return(0.0);
  208. }
  209.  
  210. #define    box_support        (0.5)
  211.  
  212. double
  213. box_filter(t)
  214. double t;
  215. {
  216.     if((t > -0.5) && (t <= 0.5)) return(1.0);
  217.     return(0.0);
  218. }
  219.  
  220. #define    triangle_support    (1.0)
  221.  
  222. double
  223. triangle_filter(t)
  224. double t;
  225. {
  226.     if(t < 0.0) t = -t;
  227.     if(t < 1.0) return(1.0 - t);
  228.     return(0.0);
  229. }
  230.  
  231. #define    bell_support        (1.5)
  232.  
  233. double
  234. bell_filter(t)        /* box (*) box (*) box */
  235. double t;
  236. {
  237.     if(t < 0) t = -t;
  238.     if(t < .5) return(.75 - (t * t));
  239.     if(t < 1.5) {
  240.         t = (t - 1.5);
  241.         return(.5 * (t * t));
  242.     }
  243.     return(0.0);
  244. }
  245.  
  246. #define    B_spline_support    (2.0)
  247.  
  248. double
  249. B_spline_filter(t)    /* box (*) box (*) box (*) box */
  250. double t;
  251. {
  252.     double tt;
  253.  
  254.     if(t < 0) t = -t;
  255.     if(t < 1) {
  256.         tt = t * t;
  257.         return((.5 * tt * t) - tt + (2.0 / 3.0));
  258.     } else if(t < 2) {
  259.         t = 2 - t;
  260.         return((1.0 / 6.0) * (t * t * t));
  261.     }
  262.     return(0.0);
  263. }
  264.  
  265. double
  266. sinc(x)
  267. double x;
  268. {
  269.     x *= M_PI;
  270.     if(x != 0) return(sin(x) / x);
  271.     return(1.0);
  272. }
  273.  
  274. #define    Lanczos3_support    (3.0)
  275.  
  276. double
  277. Lanczos3_filter(t)
  278. double t;
  279. {
  280.     if(t < 0) t = -t;
  281.     if(t < 3.0) return(sinc(t) * sinc(t/3.0));
  282.     return(0.0);
  283. }
  284.  
  285. #define    Mitchell_support    (2.0)
  286.  
  287. #define    B    (1.0 / 3.0)
  288. #define    C    (1.0 / 3.0)
  289.  
  290. double
  291. Mitchell_filter(t)
  292. double t;
  293. {
  294.     double tt;
  295.  
  296.     tt = t * t;
  297.     if(t < 0) t = -t;
  298.     if(t < 1.0) {
  299.         t = (((12.0 - 9.0 * B - 6.0 * C) * (t * tt))
  300.            + ((-18.0 + 12.0 * B + 6.0 * C) * tt)
  301.            + (6.0 - 2 * B));
  302.         return(t / 6.0);
  303.     } else if(t < 2.0) {
  304.         t = (((-1.0 * B - 6.0 * C) * (t * tt))
  305.            + ((6.0 * B + 30.0 * C) * tt)
  306.            + ((-12.0 * B - 48.0 * C) * t)
  307.            + (8.0 * B + 24 * C));
  308.         return(t / 6.0);
  309.     }
  310.     return(0.0);
  311. }
  312.  
  313. /*
  314.  *    image rescaling routine
  315.  */
  316.  
  317. typedef struct {
  318.     int    pixel;
  319.     double    weight;
  320. } CONTRIB;
  321.  
  322. typedef struct {
  323.     int    n;        /* number of contributors */
  324.     CONTRIB    *p;        /* pointer to list of contributions */
  325. } CLIST;
  326.  
  327. CLIST    *contrib;        /* array of contribution lists */
  328.  
  329. void
  330. zoom(dst, src, filterf, fwidth)
  331. Image *dst;                /* destination image structure */
  332. Image *src;                /* source image structure */
  333. double (*filterf)();            /* filter function */
  334. double fwidth;                /* filter width (support) */
  335. {
  336.     Image *tmp;            /* intermediate image */
  337.     double xscale, yscale;        /* zoom scale factors */
  338.     int i, j, k;            /* loop variables */
  339.     int n;                /* pixel number */
  340.     double center, left, right;    /* filter calculation variables */
  341.     double width, fscale, weight;    /* filter calculation variables */
  342.     Pixel *raster;            /* a row or column of pixels */
  343.  
  344.     /* create intermediate image to hold horizontal zoom */
  345.     tmp = new_image(dst->xsize, src->ysize);
  346.     xscale = (double) dst->xsize / (double) src->xsize;
  347.     yscale = (double) dst->ysize / (double) src->ysize;
  348.  
  349.     /* pre-calculate filter contributions for a row */
  350.     contrib = (CLIST *)calloc(dst->xsize, sizeof(CLIST));
  351.     if(xscale < 1.0) {
  352.         width = fwidth / xscale;
  353.         fscale = 1.0 / xscale;
  354.         for(i = 0; i < dst->xsize; ++i) {
  355.             contrib[i].n = 0;
  356.             contrib[i].p = (CONTRIB *)calloc((int) (width * 2 + 1),
  357.                     sizeof(CONTRIB));
  358.             center = (double) i / xscale;
  359.             left = ceil(center - width);
  360.             right = floor(center + width);
  361.             for(j = left; j <= right; ++j) {
  362.                 weight = center - (double) j;
  363.                 weight = (*filterf)(weight / fscale) / fscale;
  364.                 if(j < 0) {
  365.                     n = -j;
  366.                 } else if(j >= src->xsize) {
  367.                     n = (src->xsize - j) + src->xsize - 1;
  368.                 } else {
  369.                     n = j;
  370.                 }
  371.                 k = contrib[i].n++;
  372.                 contrib[i].p[k].pixel = n;
  373.                 contrib[i].p[k].weight = weight;
  374.             }
  375.         }
  376.     } else {
  377.         for(i = 0; i < dst->xsize; ++i) {
  378.             contrib[i].n = 0;
  379.             contrib[i].p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
  380.                     sizeof(CONTRIB));
  381.             center = (double) i / xscale;
  382.             left = ceil(center - fwidth);
  383.             right = floor(center + fwidth);
  384.             for(j = left; j <= right; ++j) {
  385.                 weight = center - (double) j;
  386.                 weight = (*filterf)(weight);
  387.                 if(j < 0) {
  388.                     n = -j;
  389.                 } else if(j >= src->xsize) {
  390.                     n = (src->xsize - j) + src->xsize - 1;
  391.                 } else {
  392.                     n = j;
  393.                 }
  394.                 k = contrib[i].n++;
  395.                 contrib[i].p[k].pixel = n;
  396.                 contrib[i].p[k].weight = weight;
  397.             }
  398.         }
  399.     }
  400.  
  401.     /* apply filter to zoom horizontally from src to tmp */
  402.     raster = (Pixel *)calloc(src->xsize, sizeof(Pixel));
  403.     for(k = 0; k < tmp->ysize; ++k) {
  404.         get_row(raster, src, k);
  405.         for(i = 0; i < tmp->xsize; ++i) {
  406.             weight = 0.0;
  407.             for(j = 0; j < contrib[i].n; ++j) {
  408.                 weight += raster[contrib[i].p[j].pixel]
  409.                     * contrib[i].p[j].weight;
  410.             }
  411.             put_pixel(tmp, i, k,
  412.                 (Pixel)CLAMP(weight, BLACK_PIXEL, WHITE_PIXEL));
  413.         }
  414.     }
  415.     free(raster);
  416.  
  417.     /* free the memory allocated for horizontal filter weights */
  418.     for(i = 0; i < tmp->xsize; ++i) {
  419.         free(contrib[i].p);
  420.     }
  421.     free(contrib);
  422.  
  423.     /* pre-calculate filter contributions for a column */
  424.     contrib = (CLIST *)calloc(dst->ysize, sizeof(CLIST));
  425.     if(yscale < 1.0) {
  426.         width = fwidth / yscale;
  427.         fscale = 1.0 / yscale;
  428.         for(i = 0; i < dst->ysize; ++i) {
  429.             contrib[i].n = 0;
  430.             contrib[i].p = (CONTRIB *)calloc((int) (width * 2 + 1),
  431.                     sizeof(CONTRIB));
  432.             center = (double) i / yscale;
  433.             left = ceil(center - width);
  434.             right = floor(center + width);
  435.             for(j = left; j <= right; ++j) {
  436.                 weight = center - (double) j;
  437.                 weight = (*filterf)(weight / fscale) / fscale;
  438.                 if(j < 0) {
  439.                     n = -j;
  440.                 } else if(j >= tmp->ysize) {
  441.                     n = (tmp->ysize - j) + tmp->ysize - 1;
  442.                 } else {
  443.                     n = j;
  444.                 }
  445.                 k = contrib[i].n++;
  446.                 contrib[i].p[k].pixel = n;
  447.                 contrib[i].p[k].weight = weight;
  448.             }
  449.         }
  450.     } else {
  451.         for(i = 0; i < dst->ysize; ++i) {
  452.             contrib[i].n = 0;
  453.             contrib[i].p = (CONTRIB *)calloc((int) (fwidth * 2 + 1),
  454.                     sizeof(CONTRIB));
  455.             center = (double) i / yscale;
  456.             left = ceil(center - fwidth);
  457.             right = floor(center + fwidth);
  458.             for(j = left; j <= right; ++j) {
  459.                 weight = center - (double) j;
  460.                 weight = (*filterf)(weight);
  461.                 if(j < 0) {
  462.                     n = -j;
  463.                 } else if(j >= tmp->ysize) {
  464.                     n = (tmp->ysize - j) + tmp->ysize - 1;
  465.                 } else {
  466.                     n = j;
  467.                 }
  468.                 k = contrib[i].n++;
  469.                 contrib[i].p[k].pixel = n;
  470.                 contrib[i].p[k].weight = weight;
  471.             }
  472.         }
  473.     }
  474.  
  475.     /* apply filter to zoom vertically from tmp to dst */
  476.     raster = (Pixel *)calloc(tmp->ysize, sizeof(Pixel));
  477.     for(k = 0; k < dst->xsize; ++k) {
  478.         get_column(raster, tmp, k);
  479.         for(i = 0; i < dst->ysize; ++i) {
  480.             weight = 0.0;
  481.             for(j = 0; j < contrib[i].n; ++j) {
  482.                 weight += raster[contrib[i].p[j].pixel]
  483.                     * contrib[i].p[j].weight;
  484.             }
  485.             put_pixel(dst, k, i,
  486.                 (Pixel)CLAMP(weight, BLACK_PIXEL, WHITE_PIXEL));
  487.         }
  488.     }
  489.     free(raster);
  490.  
  491.     /* free the memory allocated for vertical filter weights */
  492.     for(i = 0; i < dst->ysize; ++i) {
  493.         free(contrib[i].p);
  494.     }
  495.     free(contrib);
  496.  
  497.     free_image(tmp);
  498. }
  499.  
  500. /*
  501.  *    command line interface
  502.  */
  503.  
  504. void
  505. usage()
  506. {
  507.     fprintf(stderr, "usage: %s [-options] input.bm output.bm\n", _Program);
  508.     fprintf(stderr, "\
  509. options:\n\
  510.     -x xsize        output x size\n\
  511.     -y ysize        output y size\n\
  512.     -f filter        filter type\n\
  513. {b=box, t=triangle, q=bell, B=B-spline, h=hermite, l=Lanczos3, m=Mitchell}\n\
  514. ");
  515.     exit(1);
  516. }
  517.  
  518. void
  519. banner()
  520. {
  521.     printf("%s v%s -- %s\n", _Program, _Version, _Copyright);
  522. }
  523.  
  524. main(argc, argv)
  525. int argc;
  526. char *argv[];
  527. {
  528.     register int c;
  529.     extern int optind;
  530.     extern char *optarg;
  531.     int xsize = 0, ysize = 0;
  532.     double (*f)() = filter;
  533.     double s = filter_support;
  534.     char *dstfile, *srcfile;
  535.     Image *dst, *src;
  536.     FILE *fp;
  537.  
  538.     while((c = getopt(argc, argv, "x:y:f:V")) != EOF) {
  539.         switch(c) {
  540.         case 'x': xsize = atoi(optarg); break;
  541.         case 'y': ysize = atoi(optarg); break;
  542.         case 'f':
  543.             switch(*optarg) {
  544.             case 'b': f=box_filter; s=box_support; break;
  545.             case 't': f=triangle_filter; s=triangle_support; break;
  546.             case 'q': f=bell_filter; s=bell_support; break;
  547.             case 'B': f=B_spline_filter; s=B_spline_support; break;
  548.             case 'h': f=filter; s=filter_support; break;
  549.             case 'l': f=Lanczos3_filter; s=Lanczos3_support; break;
  550.             case 'm': f=Mitchell_filter; s=Mitchell_support; break;
  551.             default: usage();
  552.             }
  553.             break;
  554.         case 'V': banner(); exit(EXIT_SUCCESS);
  555.         case '?': usage();
  556.         default:  usage();
  557.         }
  558.     }
  559.     if((argc - optind) != 2) usage();
  560.     srcfile = argv[optind];
  561.     dstfile = argv[optind + 1];
  562.     if(((fp = fopen(srcfile, "r")) == NULL)
  563.     || ((src = load_image(fp)) == NULL)) {
  564.         fprintf(stderr, "%s: can't load source image '%s'\n",
  565.             _Program, srcfile);
  566.         exit(EXIT_FAILURE);
  567.     }
  568.     fclose(fp);
  569.     if(xsize <= 0) xsize = src->xsize;
  570.     if(ysize <= 0) ysize = src->ysize;
  571.     dst = new_image(xsize, ysize);
  572.     zoom(dst, src, f, s);
  573.     if(((fp = fopen(dstfile, "w")) == NULL)
  574.     || (save_image(fp, dst) != 0)) {
  575.         fprintf(stderr, "%s: can't save destination image '%s'\n",
  576.             _Program, dstfile);
  577.         exit(EXIT_FAILURE);
  578.     }
  579.     fclose(fp);
  580.     exit(EXIT_SUCCESS);
  581. }
  582.  
  583.