izoom.c
上传用户:weiyuanprp
上传日期:2020-05-20
资源大小:1169k
文件大小:14k
源码类别:

传真(Fax)编程

开发平台:

C/C++

  1. /* $Id: izoom.c,v 1.1.1.1 2005/11/11 21:32:03 faxguy Exp $ */
  2. /*
  3.  * Copyright (c) 1990-1996 Sam Leffler
  4.  * Copyright (c) 1991-1996 Silicon Graphics, Inc.
  5.  * HylaFAX is a trademark of Silicon Graphics
  6.  *
  7.  * Permission to use, copy, modify, distribute, and sell this software and 
  8.  * its documentation for any purpose is hereby granted without fee, provided
  9.  * that (i) the above copyright notices and this permission notice appear in
  10.  * all copies of the software and related documentation, and (ii) the names of
  11.  * Sam Leffler and Silicon Graphics may not be used in any advertising or
  12.  * publicity relating to the software without the specific, prior written
  13.  * permission of Sam Leffler and Silicon Graphics.
  14.  * 
  15.  * THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, 
  16.  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY 
  17.  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  
  18.  * 
  19.  * IN NO EVENT SHALL SAM LEFFLER OR SILICON GRAPHICS BE LIABLE FOR
  20.  * ANY SPECIAL, INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
  21.  * OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
  22.  * WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF 
  23.  * LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE 
  24.  * OF THIS SOFTWARE.
  25.  */
  26. /*
  27.  *      izoom- 
  28.  *              Magnify or minify a picture with or without filtering.  The 
  29.  *      filtered method is one pass, uses 2-d convolution, and is optimized 
  30.  *      by integer arithmetic and precomputation of filter coeffs.
  31.  *
  32.  *                              Paul Haeberli - 1988
  33.  */
  34. #include "stdio.h"
  35. #include "math.h"
  36. #include "assert.h"
  37. #include "izoom.h"
  38. #define SHIFT           12
  39. #define ONE             (1<<SHIFT)
  40. #define EPSILON         0.0001
  41. static FILTER *makefilt();
  42. static float filterrad();
  43. static float filterinteg();
  44. static float mitchell();
  45. static int (*xfiltfunc)(); 
  46. static makexmap();
  47. static xscalebuf();
  48. static addrow();
  49. static divrow();
  50. static freefilt();
  51. static applyxfilt();
  52. static int filtershape;
  53. static float blurfactor;
  54. static mults;
  55. #define GRIDTOFLOAT(pos,n)      (((pos)+0.5)/(n))
  56. #define FLOATTOGRID(pos,n)      ((pos)*(n))
  57. copyimage(getfunc,putfunc,nx,ny)
  58. int (*getfunc)(), (*putfunc)();
  59. int nx ,ny;
  60. {
  61.     int y;
  62.     short *abuf;
  63.     abuf = (short *)malloc(nx*sizeof(short));
  64.     for(y=0; y<ny; y++) {
  65.         getfunc(abuf,y);
  66.         putfunc(abuf,y);
  67.     }
  68.     free(abuf);
  69. }
  70. /*
  71.  *      general zoom follows
  72.  *
  73.  */
  74. zoom *newzoom(getfunc,anx,any,bnx,bny,filttype,blur)
  75. int (*getfunc)();
  76. int anx,any,bnx,bny;
  77. int filttype;
  78. float blur;
  79. {
  80.     zoom *z;
  81.     int i;
  82.     int xmults, ymults;
  83.     z = (zoom *)malloc(sizeof(zoom));
  84.     z->getfunc = getfunc;
  85.     z->abuf = (short *)malloc(anx*sizeof(short));
  86.     z->bbuf = (short *)malloc(bnx*sizeof(short));
  87.     z->anx = anx;
  88.     z->any = any;
  89.     z->bnx = bnx;
  90.     z->bny = bny;
  91.     z->curay = -1;
  92.     z->y = 0;
  93.     z->type = filttype;
  94.     if(filttype == IMPULSE) {
  95.         if(z->anx != z->bnx) {
  96.             z->xmap = (short **)malloc(z->bnx*sizeof(short *));
  97.             makexmap(z->abuf,z->xmap,z->anx,z->bnx);
  98.         }
  99.     } else {
  100.         filtershape = filttype;
  101.         blurfactor = blur;
  102.         if(filtershape == MITCHELL) 
  103.             z->clamp = 1;
  104.         else
  105.             z->clamp = 0;
  106.         z->tbuf = (short *)malloc(bnx*sizeof(short));
  107.         z->xfilt = makefilt(z->abuf,anx,bnx,&z->nrows);
  108.         xmults = any*mults;
  109.         z->yfilt = makefilt(0,any,bny,&z->nrows);
  110.         ymults = bnx*mults;
  111. #ifdef notdef
  112.         fprintf(stderr,"TOTAL MULTS: %dn",xmults+ymults);
  113. #endif
  114.         z->filtrows = (short **)malloc(z->nrows * sizeof(short *));
  115.         for(i=0; i<z->nrows; i++) 
  116.             z->filtrows[i] = (short *)malloc(z->bnx*sizeof(short));
  117.         z->accrow = (int *)malloc(z->bnx*sizeof(int));
  118.         z->ay = 0;
  119.     }
  120.     return z;
  121. }
  122. getzoomrow(z,buf,y)
  123. zoom *z;
  124. short *buf;
  125. int y;
  126. {
  127.     float fy;
  128.     int ay;
  129.     FILTER *f;
  130.     int i, max;
  131.     short *row;
  132.     if(y==0) {
  133.         z->curay = -1;
  134.         z->y = 0;
  135.         z->ay = 0;
  136.     }
  137.     if(z->type == IMPULSE) {
  138.         fy = GRIDTOFLOAT(z->y,z->bny);
  139.         ay = FLOATTOGRID(fy,z->any);
  140.         if(z->anx == z->bnx) {
  141.             if(z->curay != ay) {
  142.                 z->getfunc(z->abuf,ay);
  143.                 z->curay = ay;
  144.                 if(xfiltfunc) 
  145.                     xfiltfunc(z->abuf,z->bnx);
  146.             }
  147.             memcpy(buf,z->abuf,z->bnx*sizeof(short)); 
  148.         } else {
  149.             if(z->curay != ay) {
  150.                 z->getfunc(z->abuf,ay);
  151.                 xscalebuf(z->xmap,z->bbuf,z->bnx);
  152.                 z->curay = ay;
  153.                 if(xfiltfunc)
  154.                     xfiltfunc(z->bbuf,z->bnx);
  155.             }
  156.             memcpy(buf,z->bbuf,z->bnx*sizeof(short)); 
  157.         }
  158.     } else if(z->any == 1 && z->bny == 1) {
  159.             z->getfunc(z->abuf,z->ay++);
  160.             applyxfilt(z->filtrows[0],z->xfilt,z->bnx);
  161.             if(xfiltfunc)
  162.                 xfiltfunc(z->filtrows[0],z->bnx);
  163.             if(z->clamp) {
  164.                 clamprow(z->filtrows[0],z->tbuf,z->bnx);
  165.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  166.             } else {
  167.                 memcpy(buf,z->filtrows[0],z->bnx*sizeof(short)); 
  168.             }
  169.     } else {
  170.         f = z->yfilt+z->y;
  171.         max = ((int)f->dat)/sizeof(short)+(f->n-1);
  172.         while(z->ay<=max) {
  173.             z->getfunc(z->abuf,z->ay++);
  174.             row = z->filtrows[0];
  175.             for(i=0; i<(z->nrows-1); i++) 
  176.                 z->filtrows[i] = z->filtrows[i+1];
  177.             z->filtrows[z->nrows-1] = row;
  178.             applyxfilt(z->filtrows[z->nrows-1],z->xfilt,z->bnx);
  179.             if(xfiltfunc)
  180.                 xfiltfunc(z->filtrows[z->nrows-1],z->bnx);
  181.         }
  182.         if(f->n == 1) {
  183.             if(z->clamp) {
  184.                 clamprow(z->filtrows[z->nrows-1],z->tbuf,z->bnx);
  185.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  186.             } else {
  187.                 memcpy(buf,z->filtrows[z->nrows-1],z->bnx*sizeof(short)); 
  188.             }
  189.         } else {
  190.             memset(z->accrow,0,z->bnx*sizeof(int));
  191.             for(i=0; i<f->n; i++) 
  192.                 addrow(z->accrow, z->filtrows[i+(z->nrows-1)-(f->n-1)],
  193.                                                           f->w[i],z->bnx);
  194.             divrow(z->accrow,z->bbuf,f->totw,z->bnx);
  195.             if(z->clamp) {
  196.                 clamprow(z->bbuf,z->tbuf,z->bnx);
  197.                 memcpy(buf,z->tbuf,z->bnx*sizeof(short)); 
  198.             } else {
  199.                 memcpy(buf,z->bbuf,z->bnx*sizeof(short)); 
  200.             }
  201.         }
  202.     }
  203.     z->y++;
  204. }
  205. freezoom(z)
  206. zoom *z;
  207. {
  208.     int i;
  209.     if(z->type == IMPULSE) {
  210.         if(z->anx != z->bnx)
  211.             free(z->xmap);
  212.     } else {
  213.         freefilt(z->xfilt,z->bnx);
  214.         freefilt(z->yfilt,z->bny);
  215.         free(z->tbuf);
  216.         for(i=0; i<z->nrows; i++)
  217.             free(z->filtrows[i]);
  218.         free(z->filtrows);
  219.         free(z->accrow);
  220.     }
  221.     free(z->abuf);
  222.     free(z->bbuf);
  223.     free(z);
  224. }
  225. filterzoom(getfunc,putfunc,anx,any,bnx,bny,filttype,blur)
  226. int (*getfunc)(), (*putfunc)();
  227. int anx, any;
  228. int bnx, bny;
  229. int filttype;
  230. float blur;
  231. {
  232.     zoom *z;
  233.     int y;
  234.     short *buf;
  235.     buf = (short *)malloc(bnx*sizeof(short));
  236.     z = newzoom(getfunc,anx,any,bnx,bny,filttype,blur);
  237.     for(y=0; y<bny; y++) {
  238.         getzoomrow(z,buf,y);
  239.         putfunc(buf,y);
  240.     }
  241.     freezoom(z);
  242.     free(buf);
  243. }
  244. /*
  245.  *      impulse zoom utilities
  246.  *
  247.  */
  248. static makexmap(abuf,xmap,anx,bnx)
  249. short *abuf;
  250. short *xmap[];
  251. int anx, bnx;
  252. {
  253.     int x, ax;
  254.     float fx;
  255.     for(x=0; x<bnx; x++) {
  256.        fx = GRIDTOFLOAT(x,bnx);
  257.        ax = FLOATTOGRID(fx,anx);
  258.        xmap[x] = abuf+ax;
  259.     }
  260. }
  261. static xscalebuf(xmap,bbuf,bnx)
  262. short *xmap[];
  263. short *bbuf;
  264. int bnx;
  265. {
  266.     while(bnx>=8) {
  267.         bbuf[0] = *(xmap[0]);
  268.         bbuf[1] = *(xmap[1]);
  269.         bbuf[2] = *(xmap[2]);
  270.         bbuf[3] = *(xmap[3]);
  271.         bbuf[4] = *(xmap[4]);
  272.         bbuf[5] = *(xmap[5]);
  273.         bbuf[6] = *(xmap[6]);
  274.         bbuf[7] = *(xmap[7]);
  275.         bbuf += 8;
  276.         xmap += 8;
  277.         bnx -= 8;
  278.     }
  279.     while(bnx--) 
  280.         *bbuf++ = *(*xmap++);
  281. }
  282. zoomxfilt(filtfunc)
  283. int (*filtfunc)(); 
  284. {
  285.     xfiltfunc = filtfunc;
  286. }
  287. /*
  288.  *      filter zoom utilities
  289.  *
  290.  */
  291. static addrow(iptr,sptr,w,n)
  292. int *iptr;
  293. short *sptr;
  294. int w, n;
  295. {
  296.     while(n>=8) {
  297.         iptr[0] += (w*sptr[0]);
  298.         iptr[1] += (w*sptr[1]);
  299.         iptr[2] += (w*sptr[2]);
  300.         iptr[3] += (w*sptr[3]);
  301.         iptr[4] += (w*sptr[4]);
  302.         iptr[5] += (w*sptr[5]);
  303.         iptr[6] += (w*sptr[6]);
  304.         iptr[7] += (w*sptr[7]);
  305.         iptr += 8;
  306.         sptr += 8;
  307.         n -= 8;
  308.     }
  309.     while(n--) 
  310.         *iptr++ += (w * *sptr++);
  311. }
  312. static divrow(iptr,sptr,tot,n)
  313. int *iptr;
  314. short *sptr;
  315. int tot, n;
  316. {
  317.     while(n>=8) {
  318.         sptr[0] = iptr[0]/tot;
  319.         sptr[1] = iptr[1]/tot;
  320.         sptr[2] = iptr[2]/tot;
  321.         sptr[3] = iptr[3]/tot;
  322.         sptr[4] = iptr[4]/tot;
  323.         sptr[5] = iptr[5]/tot;
  324.         sptr[6] = iptr[6]/tot;
  325.         sptr[7] = iptr[7]/tot;
  326.         sptr += 8;
  327.         iptr += 8;
  328.         n -= 8;
  329.     }
  330.     while(n--)
  331.         *sptr++ = (*iptr++)/tot;
  332. }
  333. static FILTER *makefilt(abuf,anx,bnx,maxn)
  334. short *abuf;
  335. int anx, bnx;
  336. int *maxn;
  337. {
  338.     FILTER *f, *filter;
  339.     int x, n;
  340.     float bmin, bmax, bcent, brad;
  341.     float fmin, fmax, acent, arad;
  342.     int amin, amax;
  343.     float cover;
  344.     f = filter = (FILTER *)malloc(bnx*sizeof(FILTER));
  345.     *maxn = 0;
  346.     mults = 0;
  347.     for(x=0; x<bnx; x++) {
  348.         if(bnx<anx) {
  349.             brad = filterrad()/bnx;
  350.             bcent = ((float)x+0.5)/bnx;
  351.             amin = floor((bcent-brad)*anx+EPSILON);
  352.             amax = floor((bcent+brad)*anx-EPSILON);
  353.             if(amin<0)
  354.                amin = 0;
  355.             if(amax>=anx)
  356.                amax = anx-1;
  357.             f->n = 1+amax-amin;
  358.             mults += f->n;
  359.             f->dat = abuf+amin;
  360.             f->w = (short *)malloc(f->n*sizeof(short));
  361.             f->totw = 0;
  362.             for(n=0; n<f->n; n++) {
  363.                 bmin = bnx*((((float)amin+n)/anx)-bcent);
  364.                 bmax = bnx*((((float)amin+n+1)/anx)-bcent);
  365.                 cover = filterinteg(bmax,1)-filterinteg(bmin,0);  
  366.                 f->w[n] = (ONE*cover)+0.5;
  367.                 f->totw += f->w[n];
  368.             }
  369.         } else {
  370.             arad = filterrad()/anx;
  371.             bmin = ((float)x)/bnx;
  372.             bmax = ((float)x+1.0)/bnx;
  373.             amin = floor((bmin-arad)*anx+0.5+EPSILON);
  374.             amax = floor((bmax+arad)*anx-0.5-EPSILON);
  375.             if(amin<0)
  376.                amin = 0;
  377.             if(amax>=anx)
  378.                amax = anx-1;
  379.             f->n = 1+amax-amin;
  380.             mults += f->n;
  381.             f->dat = abuf+amin;
  382.             f->w = (short *)malloc(f->n*sizeof(short));
  383.             f->totw = 0;
  384.             for(n=0; n<f->n; n++) {
  385.                 acent = (amin+n+0.5)/anx;
  386.                 fmin = anx*(bmin-acent);
  387.                 fmax = anx*(bmax-acent);
  388.                 cover = filterinteg(fmax,1)-filterinteg(fmin,0);  
  389.                 f->w[n] = (ONE*cover)+0.5;
  390.                 f->totw += f->w[n];
  391.             }
  392.         }
  393.         if(f->n>*maxn)
  394.             *maxn = f->n;
  395.         f++;
  396.     }
  397.     return filter;
  398. }
  399. static freefilt(filt,n)
  400. FILTER *filt;
  401. int n;
  402. {
  403.     FILTER *f;
  404.     f = filt;
  405.     while(n--) {
  406.         free(f->w);
  407.         f++;
  408.     }
  409.     free(filt);
  410. }
  411. static applyxfilt(bbuf,xfilt,bnx)
  412. short *bbuf;
  413. FILTER *xfilt;
  414. int bnx;
  415. {
  416.     short *w;
  417.     short *dptr;
  418.     int n, val;
  419.     while(bnx--) {
  420.         if((n=xfilt->n) == 1) {
  421.             *bbuf++ = *xfilt->dat;
  422.         } else {
  423.             w = xfilt->w;
  424.             dptr = xfilt->dat;
  425.             val = 0;
  426.             n = xfilt->n;
  427.             while(n--) 
  428.                  val += *w++ * *dptr++;
  429.             *bbuf++ = val / xfilt->totw;
  430.         }
  431.         xfilt++;
  432.     }
  433. }
  434. static float filterrad()
  435. {
  436.     switch(filtershape) {
  437.         case BOX:
  438.             return 0.5*blurfactor;
  439.         case TRIANGLE:
  440.             return 1.0*blurfactor;
  441.         case QUADRATIC:
  442.             return 1.0*blurfactor;
  443.         case MITCHELL:
  444.             return 2.0*blurfactor;
  445.     }
  446.     /*NOTREACHED*/
  447. }
  448. static float quadinteg(x)
  449. float x;
  450. {
  451.    if(x<-1.0)
  452.        return 0.0;
  453.    if(x<-0.5)
  454.        return 2.0*((1.0/3.0)*(x*x*x+1.0) + x*x + x);
  455.    else
  456.        return -(2.0/3.0)*x*x*x + x + 0.5;
  457. }
  458. static float filterinteg(x,side)
  459. float x;
  460. int side;
  461. {
  462.     float val;
  463.     x = x/blurfactor;
  464.     switch(filtershape) {
  465.         case BOX:
  466.             if(x<-0.5)
  467.                 return 0.0;
  468.             else if(x>0.5)
  469.                 return 1.0;
  470.             else
  471.                 return x+0.5;
  472.         case TRIANGLE:
  473.             if(x<-1.0)
  474.                 return 0.0;
  475.             else if(x>1.0)
  476.                 return 1.0;
  477.             else if(x<0.0) {
  478.                 val = x+1.0;
  479.                 return 0.5*val*val;
  480.             } else {
  481.                 val = 1.0-x;
  482.                 return 1.0-0.5*val*val;
  483.             }
  484.         case QUADRATIC:
  485.             if(x<0.0)
  486.                 return quadinteg(x);
  487.             else
  488.                 return 1.0-quadinteg(-x);
  489.         case MITCHELL:
  490.             if(side == 0)
  491.                 return 0.0;
  492.             return mitchell(x);
  493.     }
  494.     /*NOTREACHED*/
  495. }
  496. static float p0, p2, p3, q0, q1, q2, q3;
  497. /*
  498.  * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
  499.  * SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.
  500.  */
  501. static mitchellinit(b,c)
  502. float b, c;
  503. {
  504.     p0 = (  6. -  2.*b        ) / 6.;
  505.     p2 = (-18. + 12.*b +  6.*c) / 6.;
  506.     p3 = ( 12. -  9.*b -  6.*c) / 6.;
  507.     q0 = (           8.*b + 24.*c) / 6.;
  508.     q1 = (        - 12.*b - 48.*c) / 6.;
  509.     q2 = (           6.*b + 30.*c) / 6.;
  510.     q3 = (     -     b -  6.*c) / 6.;
  511. }
  512. static float mitchell(x)        /* Mitchell & Netravali's two-param cubic */
  513. float x;
  514. {
  515.     static int firsted;
  516.     if(!firsted) {
  517.         mitchellinit(1.0/3.0,1.0/3.0);
  518.         firsted = 1;
  519.     }
  520.     if (x<-2.) return 0.0;
  521.     if (x<-1.) return 2.0*(q0-x*(q1-x*(q2-x*q3)));
  522.     if (x<0.) return 2.0*(p0+x*x*(p2-x*p3));
  523.     if (x<1.) return 2.0*(p0+x*x*(p2+x*p3));
  524.     if (x<2.) return 2.0*(q0+x*(q1+x*(q2+x*q3)));
  525.     return 0.0;
  526. }