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