home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / tools / elastic.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-06  |  19.2 KB  |  812 lines
  1. /*
  2. %  ELASTIC . C - High performance enhancement tool
  3. %
  4. %    Copyright (c) 1991    Jin, Guojun
  5. %
  6. */
  7. char    usage[]="options\n\
  8. -o option is used for PC kind computer which requires binary output    \n\
  9.     file modes. No space allowed between switch and filename.    \n\
  10. -v #    maximum value for output.    \n\
  11. -n    negative elastic operation.    \n\
  12. -b    let Elastic scaling to emphasis background. Regularly emphasis    \n\
  13.     the foreground.\n\
  14. -f #    scale factor in percentage number(in integer). default = 0%    \n\
  15.     (useful range from -255[gentle] to +450[sharp])            \n\
  16. -g[#]    base value at bottom. No output value will less than it.    \n\
  17. -a[#]    adjust float input automatically. Default=196. Smaller number    \n\
  18.     yields narrow output range.    \n\
  19. -B -F    force BYTE or FP format output.    \n\
  20. -e    emphasize in high frequency domain.    \n\
  21. -l    linear scaling    \n\
  22. -L name_of_LKT    user Look-Up Table -- 256 entries\n\
  23. -R #    interpolation Regions for both X and Y    \n\
  24. -Rx #    X interpolation regions.\n\
  25. -Ry #    Y interpolation regions.\n\
  26. -S    treat all slices as Single 3D Image (one frame).\n\
  27. -t [#n [#v]]    \n\
  28.     clip top #n(default=1) value to 91%(default) or given value #v.    \n\
  29. -z    option to count zero valued pixels as frequency factor(for -e).    \n\
  30. -r    option only for FAST_Version to get relaxed elastic scale.    \n\
  31. -E    modify Edges on right and bottom\n\
  32. -I    Interpolation in main peak area for BYTE format output. [2 x 2]    \n\
  33. -M    Message switch on.    \n\
  34. -T    generate Table for plot    \n\
  35. [<] in_file [> [-o] output_file]\n", *u_lkt;
  36.  
  37. /* Note:
  38. %    input can be any format.
  39. %    FP input can only output either BYTE or FP without interpolation.
  40. %
  41. % compile: cc -O -o elastic elastic.c -lscs3 -lccs -lhips -lrle -ltiff -lm
  42. %
  43. % AUTHOR:    Jin Guojun - LBL    1/3/91
  44. */
  45.  
  46. #include "header.def"
  47. #include "imagedef.h"
  48. #include <math.h>    /* new_curve required    */
  49.  
  50. #define    ADJ_fact    196
  51. #define    GValue()    arget(argc, argv, &i, &j)
  52.  
  53. /*    GLOBAL variables    */
  54.  
  55. bool    adj, bgrd, flr, interpolate, neg, Msg, top, Table;
  56. #ifdef    FAST
  57. bool    revs=0;
  58. #define    VER    "Fast"        /* strange function    */
  59. #else
  60. #define    VER    "Regular"    /* spring law    */
  61. #endif
  62. #ifndef    VSC_BASE
  63. #define    VSC_BASE    .01
  64. #endif
  65.  
  66. bool    D3, edge, hfeq, zflag, zcount;
  67. int    regions, rgn_w, rgn_h, xfrac, yfrac;
  68. float    h_avg, hint, maxout, rel_val;
  69. MType    fsize, sc=0, topv, maxval;
  70. U_IMAGE    uimg;
  71.  
  72. #define    cols    uimg.width
  73. #define    rows    uimg.height
  74. #define    nfrm    uimg.frames
  75. #define    ibuf    uimg.src
  76. #define    obuf    uimg.dest
  77.  
  78.  
  79. main(argc, argv)
  80. int    argc;
  81. char **    argv;
  82. {
  83. #define    p    i
  84.  
  85. byte    *obp;
  86. int    *ip, *tmpip, ffac;
  87. float    *ofp, fmin, fmin_1, *fip;
  88. InterpMap    *IM;    /* Interpolation Matrix */
  89. MType    frmp, i, j, x_regions=2, y_regions=2;
  90. LKT*    lktp;
  91.  
  92. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, uimg.o_form = -1, *argv, "S20-1");
  93.  
  94. for (i=1; i<argc; i++)
  95.     if (*argv[i] == '-')    {
  96.     j = 1;
  97.     switch(argv[i][j++])
  98.     {
  99.     case 'a':if ((adj=GValue()) == 0)
  100.             adj=ADJ_fact;    break;
  101.     case 'b':    bgrd++;    break;
  102.     case 'e':    bgrd--;    break;
  103.     case 'f':    sc = GValue();    break;
  104.     case 'g':    flr = GValue();
  105.         if (!flr)    flr=10;    break;
  106.     case 'E':    edge++;    break;
  107.     case 'L':
  108.         if (avset(argc, argv, &i, &j, 1))
  109.             u_lkt = argv[i] + j;
  110.         break;
  111.     case 'R':{
  112.     register int    who=argv[i][j], num=GValue();
  113.         switch(who){
  114.         case 'x':    x_regions = num;    break;
  115.         case 'y':    y_regions = num;    break;
  116.         default:    x_regions = y_regions = num;
  117.         }
  118.     }
  119.     case 'I':    interpolate++;    break;
  120.     case 'B':    uimg.o_form=IFMT_BYTE;    break;
  121.     case 'F':    uimg.o_form=IFMT_FLOAT;    break;
  122.     case 'M':    message("%s: %s\n", *argv, VER);
  123.         Msg++;    break;
  124. #ifdef    _DEBUG_
  125.     case 'D':    debug++;    break;
  126. #endif
  127.     case 'l':    bgrd=2;    break;    /* linear op    */
  128.     case 'n':    neg++;    break;
  129. #ifdef    FAST
  130.     case 'r':    revs++;    break;
  131. #endif
  132.     case 'S':    D3++;    break;
  133.     case 't':    top = GValue()-1;
  134.             topv = GValue();
  135.         break;
  136.     case 'T':    Table++;    break;
  137.     case 'v':    maxval = GValue();    break;
  138.     case 'z':
  139.     case 'Z':    zflag++;    break;
  140.     case 'o':
  141.         if (avset(argc, argv, &i, &j, 1) &&
  142.             freopen(argv[i]+j, "wb", stdout))    break;
  143.         message("output file - %s", argv[i]);
  144.     default:
  145. info:        usage_n_options(usage, i, argv[i]);
  146.     }
  147.     }else if ((in_fp=freopen(argv[i], "rb", stdin)) == 0)
  148.         syserr("input file - %s", argv[i]);
  149.  
  150. io_test(fileno(in_fp), goto    info);
  151.  
  152. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  153. fsize = cols * rows;
  154.  
  155. if (D3)    fsize *= nfrm,    nfrm = 1;
  156.  
  157. if (uimg.o_form==IFMT_BYTE)
  158.     uimg.pxl_out = 1;
  159. else if (uimg.o_form==IFMT_FLOAT)
  160.     uimg.pxl_out = sizeof(float);
  161. else    uimg.o_form = uimg.in_form,
  162.     uimg.pxl_out = uimg.pxl_in;
  163.  
  164. if (!maxval)
  165.     switch(uimg.pxl_in){
  166.     case 1:    maxval=255;    break;
  167.     default:maxval=65535;    break;
  168.     }
  169. if (top<0 || top>maxval)    top = 1;
  170. if (top){
  171.     top = maxval - top;
  172.     if (!topv || topv>=maxval)
  173.         topv = maxval * .91;
  174. }
  175. if (bgrd<0 && uimg.in_form>IFMT_BYTE)
  176.     if (uimg.in_form==IFMT_FLOAT)    bgrd = 1;
  177.     else    interpolate++;
  178. interpolate &= (uimg.o_form<=IFMT_LONG && !u_lkt);
  179.  
  180. ibuf = nzalloc(fsize, uimg.pxl_in, "ibuf");
  181. uimg.hist = (int*)nzalloc(maxval+1,sizeof(*(uimg.hist)), "hist");
  182.  
  183. if (!interpolate)
  184.     lktp = (LKT*)(uimg.lkt = zalloc(maxval, sizeof(LKT), "lkt"));
  185. if (uimg.in_form != IFMT_LONG)
  186.     tmpip = (int*)zalloc(fsize, sizeof(*tmpip), "tmp");
  187. if (u_lkt){
  188. register int    t_l;
  189. FILE    *lkt_fp = fopen(u_lkt, "r");
  190.     if (!lkt_fp)    syserr("user lkt %s", u_lkt);
  191.     t_l = getw(lkt_fp);
  192.     if (t_l > maxval)    t_l = maxval;
  193.     for (i=0; i<t_l; i++)
  194.         fscanf(lkt_fp, "%f", &lktp[i]);
  195. }
  196.  
  197. if (uimg.pxl_out==1)
  198.     maxout = 255;
  199. else    maxout = maxval;
  200.  
  201. if (!Table)    (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  202.  
  203. if (interpolate){
  204.     if (y_regions<2 || y_regions>(rows>>4) || y_regions>16)
  205.         y_regions = 4;
  206.     if (x_regions<2 || x_regions>(cols>>4) || x_regions>16)
  207.         x_regions = 4;
  208.     message("x regions=%d,    y regions=%d\n", x_regions, y_regions);
  209.     rgn_h = rows / y_regions;
  210.     rgn_w = cols / x_regions;
  211.     xfrac = cols - rgn_w * x_regions;
  212.     yfrac = rows - rgn_h * y_regions;
  213.     regions = (x_regions + (xfrac!=0)) * (y_regions + (yfrac!=0));
  214.     IM = (InterpMap*)nzalloc(regions, sizeof(*IM), "IM");
  215.     for (i=0; i<regions; i++)
  216.         IM[i].lkt = (LKT*)zalloc(maxval, sizeof(LKT), "lkts");
  217. }
  218.  
  219. if (bgrd>=0 && !interpolate)
  220.     obuf = nzalloc(fsize, uimg.pxl_out, "obuf");
  221. else    obuf = ibuf;
  222.  
  223. for (frmp=0; frmp < nfrm; frmp++)
  224. {
  225. int    min=uimg.o_form,    max=uimg.pxl_out;
  226.  
  227.     if (uimg.in_type==FITS)
  228.         uimg.o_form = uimg.in_form,    uimg.pxl_out = uimg.pxl_in;
  229.     (*uimg.std_swif)(FI_LOAD_FILE, &uimg, NULL, No);
  230.     uimg.o_form = min,    uimg.pxl_out = max;
  231.  
  232.     if (bgrd>=0 || interpolate){
  233.     ip = tmpip;
  234.     switch(uimg.in_form){
  235.     case IFMT_BYTE:    obp = (byte*)ibuf;
  236.         for (i=fsize; i--;)    *ip++ = *obp++;
  237.         break;
  238.     case IFMT_SHORT: {
  239.     register short    *osp = (short*)ibuf;
  240.         for (i=fsize; i--;)    *ip++ = *osp++;
  241.     }    break;
  242.     case IFMT_LONG:    tmpip = (int*)ibuf;
  243.         break;
  244.     case IFMT_FLOAT:
  245.         maxval = ffind_min_max(ibuf, &fmin, &fmin_1, &max, &ffac);
  246.         uimg.lkt = zalloc(maxval, (MType)sizeof(LKT), "fp_lkt");
  247.         if (!u_lkt)    new_curve(uimg.lkt, maxval, 0, max, fsize);
  248.         fip = (float*)ibuf;
  249.         switch(uimg.o_form){
  250.         case IFMT_BYTE:
  251.             for (p=fsize, obp=(byte*)obuf; p--; fip++)
  252.             *obp++ = lktp[(int)(ffac*(*fip-fmin) * fmin_1)];
  253.         break;
  254.         case IFMT_FLOAT:
  255.             for (p=fsize, ofp=(float*)obuf; p--; fip++)
  256.             *ofp++ = lktp[(int)(ffac*(*fip-fmin) * fmin_1)];
  257.         }
  258.         free(uimg.lkt);
  259.         goto    wout;
  260.     default:syserr("Only accept B S I F format");
  261.     }
  262.     ip = tmpip;
  263.     }
  264.     else    ip = (int*)ibuf;    /* histo equalization */
  265.     if (interpolate)
  266.     interpolation(ip, obuf, x_regions, y_regions, IM);
  267.     else{
  268.     maxval = Find_min_max(ip, &min, &max);
  269.  
  270.     /* building new histogram table    */
  271.     if (!u_lkt)    new_curve(uimg.lkt, maxval, min, max, fsize);
  272.  
  273.     /* generate new histogram    */
  274.     if (bgrd<0){
  275.     register byte    *bp = PrtCAST ip;
  276.         for (p=fsize; p--;)
  277.             *bp = lktp[*bp++];
  278.     }
  279.     else switch (uimg.o_form)
  280.     {
  281.     case IFMT_BYTE:
  282.         for (p=fsize, obp=(byte*)obuf; p--;)
  283.             *obp++ = lktp[*ip++ - min] + flr;
  284.         break;
  285.     case IFMT_SHORT: {
  286.     register short    *osp = (short*)obuf;
  287.         for (p=fsize; p--;)
  288.             *osp++ = lktp[*ip++ - min] + flr;
  289.         } break;
  290.     case IFMT_LONG:{
  291.     register int    *iobp=(int*)obuf;
  292.         for (p=fsize; p--;)
  293.             *iobp++ = lktp[*ip++ - min] + flr;
  294.         } break;
  295.     case IFMT_FLOAT:
  296.     default:for (p=fsize, ofp=(float*)obuf; p--;)
  297.             *ofp++ = lktp[*ip++ - min] + flr;
  298.     }
  299.     }
  300. wout:    i = (*uimg.std_swif)(FI_SAVE_FILE, &uimg, uimg.dest, uimg.save_all=0);
  301. }
  302. return    (D3);
  303. }
  304.  
  305.  
  306. /* calculate the histogram */
  307. histogram(vbuf, histp, hsize, bytes)
  308. void    *vbuf;
  309. register int    *histp;
  310. {
  311. register int    i;
  312.  
  313. for (i=maxval; i--;)    histp[i]=0;
  314. switch(bytes){
  315. case 1:    {
  316. register byte    *bp = vbuf;
  317.     for (p=hsize; p--;)    ++histp[*bp++];
  318.     }break;
  319. case 2:    {
  320. register short    *sp = vbuf;
  321.     for (p=hsize; p--;)    ++histp[*sp++];
  322.     }break;
  323. case 4:    {
  324. register int    *ip = vbuf;
  325.     for (p=hsize; p--;)    ++histp[*ip++];
  326.     }break;
  327. default:    syserr("no histo for this format %d", uimg.in_form);
  328. }
  329. }
  330.  
  331. /*==============================================*
  332. *    computing the min & max for integer    *
  333. *    & trans float to integer for scaling    *
  334. *==============================================*/
  335. find_min_max(bufp, imp, r_width, r_height)
  336. register int    *bufp;
  337. InterpMap    *imp;
  338. {
  339. register unsigned    i, j=0;
  340. register int    *histp=uimg.hist;
  341. for (i=maxval; i--;)
  342.     histp[i] = j;
  343. for (i=r_height; i--; bufp+=cols)
  344.     for (j=r_width; j--;)
  345.     histp[bufp[j]]++;
  346. for (i=0, j=maxval; histp[i]==0 && i<j; i++);
  347. while (histp[--j]==0);
  348. imp->min = i;
  349. imp->max = j;
  350. return    imp->diff = j - i + 1;
  351. }
  352.  
  353. Find_min_max(buf, min, max)
  354. void    *buf;
  355. int    *min, *max;
  356. {
  357. register int    i, j;
  358.  
  359. if (uimg.pxl_in==4){
  360. register unsigned    min=65536, max=0;
  361. register int    *ip = buf;
  362.     for (i=fsize; i--;){
  363.         j = *ip++;
  364.         if (j>max)    max = j;
  365.         else if (j<min)    min = j;
  366.     }
  367.     i = min,    j = max;
  368. }
  369. else {
  370. register int    *histp=uimg.hist;
  371.     histogram(buf, histp, fsize, bgrd<0?uimg.pxl_in:sizeof(MType));
  372.     for (i=0, j=maxval; histp[i]==0 && i<j; i++);
  373.     while (histp[--j]==0);
  374. }
  375. *min = i;
  376. *max = j;
  377. return    j - i + 1;
  378. }
  379.  
  380. ffind_min_max(bufp, fmin, fmin_1, max, ffac)
  381. register float    *bufp, *fmin, *fmin_1;
  382. register int    *max, *ffac;
  383. {
  384. register int    i;
  385. register float    fmax = -3.4e38;
  386. *fmin = 3.4e38;
  387.  
  388. for (i=0; i<fsize; i++, bufp++)
  389.     if (*bufp > fmax)    fmax = *bufp;
  390.     else if (*bufp < *fmin)    *fmin = *bufp;
  391. if (fabs(fmax) > 1 && fabs(*fmin) > 1)
  392.     *fmin_1 = 1;
  393. else    if (!*fmin)    *fmin_1 = 1e7;
  394.     else    *fmin_1 = fabs(1. / *fmin);
  395. fmax = fabs(*fmin_1 * (fmax - *fmin));
  396.  
  397. if (adj)    *ffac = adj / fmax;
  398. else if (fmax < 64)    *ffac = ADJ_fact /fmax;
  399.     else    *ffac = 1;
  400. #ifdef    _DEBUG_
  401. message("fmin=%f, fm_1=%f, fmax=%f, fac=%d\n", *fmin, *fmin_1, fmax, *ffac);
  402. #endif
  403. return    *max = *ffac * fmax;
  404. }
  405.  
  406. /*======================================================*
  407. *    computing the average new value of pixels    *
  408. *======================================================*/
  409.  
  410. #ifdef    Fast
  411. new_curve(lkt, maxdiff, min, max)
  412. LKT*    lkt;
  413. int    maxdiff;
  414. register int    min, max;
  415. {
  416. int    i;
  417. register float    scale, vsf;
  418.  
  419. #ifdef    _DEBUG_
  420. if (Msg)
  421.     message("in -> min=%d, max=%d, maxdiff=%d : max_out=%f",
  422.     min, max, maxdiff, maxout);
  423. #endif
  424.     vsf = 0.01 * sc * max;
  425.     scale = maxout / (vsf + ((min + max)*maxdiff >> 1));
  426. #ifdef    _DEBUG_
  427.     if(Msg)    message("    scale=%f, sc=%f\n", scale, vsf);
  428. #endif
  429.     if (revs)
  430.        if (!bgrd)
  431.         for(i=rel_val=0; i<maxdiff; i++)
  432.         {
  433.         rel_val += (maxdiff-i+vsf) * scale;
  434.         if (rel_val > hint)    lkt[i]=hint;
  435.         else    lkt[i] = rel_val;
  436.         }
  437.        else for (i=maxdiff, rel_val=hint; i; i--)
  438.         {
  439.         lkt[i] = rel_val;
  440.         rel_val -= (i+vsf)*scale;
  441.         if (rel_val<0 && !uimg.o_form)    rel_val=0;
  442.         }
  443.     else if (neg)
  444.        if (!bgrd)
  445.         for(i=maxdiff, rel_val=0; i; i--)
  446.         {
  447.         rel_val += (maxdiff-i+vsf) * scale;
  448.         if (rel_val > maxout)    lkt[i]=maxout;
  449.         else    lkt[i] = rel_val;
  450.         }
  451.        else for (i=0, rel_val=maxout; i<maxdiff; i++)
  452.         {
  453.         lkt[i] = rel_val;
  454.         rel_val -= (i+vsf)*scale;
  455.         if (rel_val<0 && !uimg.o_form)    rel_val=0;
  456.         }
  457.     else if (!bgrd)
  458.         for(i=rel_val=0; i<maxdiff; i++)
  459.         {
  460.         rel_val += (i+vsf) * scale;
  461.         if (rel_val > maxout)    lkt[i]=maxout;
  462.         else    lkt[i] = rel_val;
  463.         }
  464.        else for (i=maxdiff, rel_val=maxout; i; i--)
  465.         {
  466.         lkt[i] = rel_val;
  467.         rel_val -= (maxdiff-i+vsf)*scale;
  468.         if (rel_val<0 && !uimg.o_form)    rel_val=0;
  469.         }
  470. #else
  471. /*===============================
  472. *    maxout is float point    *
  473. *    maxdiff is max-min+1    *
  474. *    the 1 is a cell for max    *
  475. ===============================*/
  476. new_curve(lkt, maxdiff, min, max, hsize)
  477. LKT*    lkt;
  478. MType    maxdiff;
  479. register int    min, max;
  480. {
  481. int    i;
  482.  
  483. if (bgrd < 0){
  484. int    z, incr;
  485. register int    *histp = uimg.hist;
  486.  
  487.     /* calculate the average number of pixels per bin */
  488.         if (!zflag){    zcount = histp[0];    histp[0] = 0;    }
  489.         h_avg = (float) (hsize - zcount)/maxout;
  490.  
  491.     /*==============================================================*
  492.     *    Re-calculate the average number of pixels per bin    *
  493.     *==============================================================*/
  494.  
  495.         for (z = rel_val = hint = 0; z < maxdiff; z++) {
  496.             hint += histp[z];
  497.             incr = hint/h_avg;
  498.             lkt[z] = rel_val + (incr>>2);
  499.             rel_val += incr;
  500.             hint -= incr*h_avg;
  501. #ifdef    _DEBUG_
  502.         if(Msg)    message("%6.2f(%4d)    ", lkt[z], histp[z]);
  503. #endif
  504.         };
  505.         if (Table)    GTable(maxval, uimg.lkt);
  506. }
  507. else if (bgrd==2){
  508. register float    scale = (maxout - flr) / (max - min);
  509.     for (i=0; i < maxdiff; i++)
  510.         lkt[i] = i * scale + flr;
  511. }
  512. else{
  513. register float    scale, vsc=VSC_BASE*sc/max + 1;
  514.     if (vsc<=0.)    vsc = VSC_BASE;
  515.     {
  516.     register double    tmp;
  517.     if (vsc != 1.)
  518.     {
  519.         for (i=tmp=maxdiff; i--;)    tmp = tmp*vsc + i;
  520.         scale = (maxout-flr) / tmp;
  521.     }
  522.     else    scale = (maxout-flr) / ((maxdiff-1)*maxdiff >> 1);
  523.     if (Msg)
  524.     message("min=%d, max=%d, flr=%d, top=%.3f\n", min, max, flr, maxout);
  525. #    ifdef    _DEBUG_
  526.     message("Scale=%f, C=%f, sc=%f\n", scale, tmp, vsc);
  527. #    endif
  528.     }
  529.     if (neg)
  530.        if (!bgrd)
  531.         for(i=maxdiff, rel_val=0; i--;)
  532.         {
  533.         rel_val += (maxdiff-i) * (scale*=vsc);
  534.         if (rel_val > maxout)    lkt[i]=maxout;
  535.         else    lkt[i] = rel_val;
  536.         }
  537.        else for (i=0, rel_val=maxout; i<maxdiff; i++)
  538.         {
  539.         lkt[i] = rel_val;
  540.         rel_val -= i * (scale*=vsc);
  541.         if (rel_val<0 && !uimg.o_form)    rel_val=0;
  542.         }
  543.     else if (!bgrd)
  544.         for(i=1, rel_val=0; i<maxdiff; i++)
  545.         {
  546.         rel_val += i * (scale*=vsc);
  547.         lkt[i] = rel_val;
  548.         }
  549.        else for (i=maxdiff, rel_val=maxout; i--;)
  550.         {
  551.         lkt[i] = rel_val;
  552.         rel_val -= (maxdiff-i) * (scale*=vsc);
  553.         }
  554. #endif    Fast
  555.     i = maxdiff - 1;
  556.     if (uimg.in_form==IFMT_BYTE){
  557.         if (*lkt < 0){
  558. #    ifdef    _DEBUG_
  559.             msg("lkt[0]=%d\n", *lkt);
  560. #    endif
  561.             *lkt=0;
  562.         }
  563.         else if (*lkt > 255){
  564. #    ifdef    _DEBUG_
  565.             msg("lkt[0]=%d\n", *lkt);
  566. #    endif
  567.             *lkt=255;
  568.         }
  569.         if (lkt[i] > 255)    lkt[i]=255;
  570.         else if (lkt[i]<0)    lkt[i]=0;
  571.     }
  572. }
  573. if (top)
  574.     for (;lkt[i]>top; i--)
  575.     lkt[i] = topv;
  576. if (Table)    GTable(maxdiff, lkt);
  577. if (Msg)    dump_lkt(maxdiff, lkt);
  578. }
  579.  
  580. dump_lkt(n, lkt)
  581. LKT*    lkt;
  582. {
  583. register int    i;
  584. for (i=0; i<n;){
  585.     message("%10.3f", lkt[i]);    if (!(++i & 7))    mesg("\n");
  586. }
  587. msg("\n%d items\n", i);
  588. }
  589.  
  590. GTable(max, lkt)
  591. LKT*    lkt;
  592. {
  593. register int i;
  594. register float *lp=lkt;
  595. for (i=0; i<max; i++)
  596.     printf("%d %f\n", i, *lp++);
  597. exit(0);
  598. }
  599.  
  600. interpolation(in, out, x_regions, y_regions, imp)
  601. int    *in;
  602. byte    *out;
  603. InterpMap    *imp;
  604. {
  605. int    i, j, maxdiff, x, y, mx, my, region;
  606. float    xinc, yinc, xfracinc, yfracinc, xrate, yrate, xcomp, ycomp;
  607. InterpMap    *mlu, *mld, *mru, *mrd, *imu, *imd;
  608.  
  609. {
  610. register int    *tmpp=in, diff;
  611. for (i=region=maxdiff=0; i<y_regions; i++, tmpp+=(rgn_h-1)*cols){
  612.     for (j=0; j<x_regions; j++, tmpp+=rgn_w, region++){
  613.     diff=find_min_max(tmpp, imp+region, rgn_w, rgn_h);
  614.     if (maxdiff < diff)
  615.         maxdiff = diff;
  616.     /* building new histogram table    */
  617.     new_curve(imp[region].lkt, imp[region].diff,
  618.         imp[region].min, imp[region].max, rgn_h*rgn_w);
  619.     }
  620.     if (xfrac){
  621.     diff=find_min_max(tmpp, imp+region, xfrac, rgn_h);
  622.     if (maxdiff < diff)
  623.         maxdiff = diff;
  624.     new_curve(imp[region].lkt, imp[region].diff,
  625.         imp[region].min, imp[region].max, rgn_h*xfrac);
  626.     region++;
  627.     tmpp += xfrac;
  628.     }
  629. }
  630. if (yfrac){
  631.     for (j=0; j<x_regions; j++, region++, tmpp+=rgn_w){
  632.     diff=find_min_max(tmpp, imp+region, rgn_w, yfrac);
  633.     if (maxdiff < diff)
  634.         maxdiff = diff;
  635.     new_curve(imp[region].lkt, imp[region].diff,
  636.         imp[region].min, imp[region].max, yfrac*rgn_w);
  637.     }
  638.     if (xfrac){
  639.     diff=find_min_max(tmpp, imp+region, xfrac, yfrac);
  640.     if (maxdiff < diff)
  641.         maxdiff = diff;
  642.     new_curve(imp[region].lkt, imp[region].diff,
  643.         imp[region].min, imp[region].max, yfrac*xfrac);
  644.     tmpp += xfrac;
  645.     }
  646. }
  647. }
  648.  
  649. {
  650. register int    xd=x_regions + (xfrac!=0);
  651.  
  652. for (i=regions-(yfrac?xd:0); i--;)    /* difference interpolation */
  653.     if (imp[i].diff<maxdiff && (!xfrac || (i+1)%xd)){
  654. #ifdef    _DEBUG_
  655.     message("lkt%d = %d\n", i, imp[i].diff);
  656. #endif
  657.     for (j=imp[i].diff; j<maxdiff; j++)
  658.         imp[i].lkt[j] = imp[i].lkt[j-1];
  659.     }
  660. }
  661. /* generate new histogram    */
  662. xinc = 1. / rgn_w;
  663. yinc = 1. / rgn_h;
  664. xfracinc = 1. / xfrac;
  665. yfracinc = 1. / yfrac;
  666. mx = x_regions + (xfrac != 0);
  667. for (i=0; i<y_regions; i++){
  668.     yrate = 1.;
  669.     my = i * mx;
  670.     if (i+1==y_regions)
  671.     mx = -mx;
  672.     imu = imp + my;
  673.     imd = imu + mx;
  674.     for (y=0; y<rgn_h; y++, yrate-=yinc){
  675.     ycomp = 1. - yrate;
  676.     for (j=0; j<x_regions; j++){
  677.         xrate = 1.;
  678.         mlu = imu + j;
  679.         mld = imd + j;
  680.         if (j+1==x_regions){
  681.         mru = mlu - 1;
  682.         mrd = mld - 1;
  683.         }
  684.         else{
  685.         mru = mlu + 1;
  686.         mrd = mld + 1;
  687.         }
  688.         for (x=0; x<rgn_w; x++, xrate-=xinc){
  689.         register int    A, B, C, D, where = *in++;
  690.         A = where - mlu->min;
  691.         B = where - mld->min;
  692.         C = where - mru->min;
  693.         D = where - mrd->min;
  694.         if (A < 0)    A = 0;
  695.         if (B < 0)    B = 0;
  696.         if (C < 0)    C = 0;
  697.         if (D < 0)    D = 0;
  698.         xcomp = 1. - xrate;
  699.         *out++ = xrate * (mlu->lkt[A] * yrate + mld->lkt[B] * ycomp) +
  700.             xcomp * (mru->lkt[C] * yrate + mrd->lkt[D] * ycomp);
  701.         }
  702.     }
  703.     if (xfrac){
  704. /*
  705.         xrate = 1.;
  706.         mlu = imu + j;
  707.         mld = imd + j;
  708.         mru = mlu - 1;
  709.         mrd = mld - 1;
  710.         for (x=0; x<xfrac; x++, xrate-=xfracinc){
  711.         register int    A, B, C, D, where = *in++;
  712.         A = where - mlu->min;
  713.         B = where - mld->min;
  714.         C = where - mru->min;
  715.         D = where - mrd->min;
  716.         if (A < 0)    A = 0;
  717.         if (B < 0)    B = 0;
  718.         if (C < 0)    C = 0;
  719.         if (D < 0)    D = 0;
  720.         xcomp = 1. - xrate;
  721.         *out++ = xrate * (mlu->lkt[A] * yrate + mld->lkt[B] * ycomp) +
  722.             xcomp * (mru->lkt[C] * yrate + mrd->lkt[D] * ycomp);
  723.         }
  724. */
  725.         for (x=xfrac; x--; *out++ = *in++);
  726.     }
  727.     }
  728. }
  729. if (edge)
  730. if (yfrac){
  731.     yrate = 1.;
  732.     mx = x_regions + (xfrac != 0);
  733.     my = i * mx;
  734.     imu = imp + my;
  735.     imd = imu - mx;
  736.     for (y=0; y<yfrac; y++, yrate-=yfracinc){
  737.     ycomp = 1. - yrate;
  738.     for (j=0; j<x_regions; j++){
  739.         xrate = 1.;
  740.         mlu = imu + j;
  741.         mld = imd + j;
  742.         if (j+1==x_regions && !xfrac){
  743.         mru = imu + j-1;
  744.         mrd = imd + j-1;
  745.         }
  746.         else{
  747.         mru = imu + j+1;
  748.         mrd = imd + j+1;
  749.         }
  750.         for (x=0; x<rgn_w; x++, xrate-=xinc){
  751.         register int    A, B, C, D, where = *in++;
  752.         A = where - mlu->min;
  753.         B = where - mld->min;
  754.         C = where - mru->min;
  755.         D = where - mrd->min;
  756.         if (A < 0)    A = 0;
  757.         if (B < 0)    B = 0;
  758.         if (C < 0)    C = 0;
  759.         if (D < 0)    D = 0;
  760.         xcomp = 1. - xrate;
  761.         *out++ = xrate * (mlu->lkt[A] * yrate + mld->lkt[B] * ycomp) +
  762.             xcomp * (mru->lkt[C] * yrate + mrd->lkt[D] * ycomp);
  763.         }
  764.     }
  765.     if (xfrac){
  766.         xrate = 1.;
  767.         mlu = imu + j;
  768.         mld = imd + j;
  769.         mru = imu + j-1;
  770.         mrd = imd + j-1;
  771.         for (x=0; x<xfrac; x++, xrate-=xfracinc){
  772.         register int    A, B, C, D, where = *in++;
  773.         A = where - mlu->min;
  774.         B = where - mld->min;
  775.         C = where - mru->min;
  776.         D = where - mrd->min;
  777.         if (A < 0)    A = 0;
  778.         if (B < 0)    B = 0;
  779.         if (C < 0)    C = 0;
  780.         if (D < 0)    D = 0;
  781.         xcomp = 1. - xrate;
  782.         *out++ = xrate * (mlu->lkt[A] * yrate + mld->lkt[B] * ycomp) +
  783.             xcomp * (mru->lkt[C] * yrate + mrd->lkt[D] * ycomp);
  784.         }
  785.     }
  786.     }
  787. }
  788. else {
  789. /*    if (xfrac){
  790.         xrate = 1.;
  791.         mlu = mru = imu + j-1;
  792.         mld = mrd = imd + j-1;
  793.         for (x=0; x<xfrac; x++, xrate-=xfracinc){
  794.         register int    where = *in++;
  795.         xcomp = 1. - xrate;
  796.         *out++ = xrate * (mlu->lkt[where - mlu->min] * yrate +
  797.                 mld->lkt[where - mld->min] * ycomp) +
  798.             xcomp * (mru->lkt[where - mru->min] * yrate +
  799.                 mrd->lkt[where - mrd->min] * ycomp);
  800.         }
  801.     }
  802.     }
  803. }
  804. */
  805. if (yfrac)
  806.     for (j=cols*yfrac; i--;)
  807.         *out++ = *in++;
  808.  
  809.  
  810. }
  811. }
  812.