home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / vfft / dcvtob.c next >
Encoding:
C/C++ Source or Header  |  1992-06-01  |  11.5 KB  |  426 lines
  1. /*
  2. %    DCVtoB . C
  3. %
  4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  5.  
  6. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  7. Permission is granted to reproduce this software for non-commercial
  8. purposes provided that this notice is left intact.
  9.  
  10. It is acknowledged that the U.S. Government has rights to this software
  11. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  12. and the University of California.
  13.  
  14. This software is provided as a professional and academic contribution
  15. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  16. with no warranties of any kind whatsoever, no support, no promise of
  17. updates, or printed documentation. By using this software, you
  18. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  19. University of California shall have no liability with respect to the
  20. infringement of other copyrights by any part of this software.
  21.  
  22. For further information about this notice, contact William Johnston,
  23. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  24. (wejohnston@lbl.gov)
  25.  
  26. For further information about this software, contact:
  27.     Jin Guojun
  28.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  29.     g_jin@lbl.gov
  30.  
  31. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  32. %
  33. %    This program convert following format image to Byte image:
  34. %
  35. %        float complex
  36. %        double, double complex
  37. %        float vfft 2D, float vfft 3D
  38. %        double vfft 2D, double vfft 3D
  39. %
  40. %    Also    converting VFFT and double VFFT to corresponding complex FFT.
  41. */
  42. char    usage[]="options\n\
  43. [-div #]    maximum input value divident for auto scale.    \n\
  44.         The large value results a bright picture    \n\
  45. [-fflip]    do fast quad flip on regular complex format for view only\n\
  46. [-flip ]    do regular quadrant flip on regular complex format    \n\
  47. [-maxi #]    specify the maximum input value(range)    \n\
  48. [-Complex]    VFFT to Complex (FFT)        \n\
  49. [-Float]    output Float format instead of Bytes    \n\
  50. [-nonscale]    output non scaled conversion    \n\
  51. [-regular]    output original VFFT spectrum, not quad flipped    \n\
  52. [-selfscale]    re-scaling for each frame    \n\
  53.         (regular scale calculation done at first frame)    \n\
  54. [<] input [> byte_image]\n";
  55. /*
  56. % compile:
  57. %    cc -O -o DEST/dcvtob dcvtob.c -lscs3 -lccs -lhips -lrle -ltiff -lm
  58. %
  59. % AUTHOR:    Guojun Jin - LBL    5/1/91
  60. */
  61.  
  62. #include <math.h>
  63. #include "complex.h"
  64. #include "header.def"
  65. #include "imagedef.h"
  66.  
  67. U_IMAGE    uimg;
  68.  
  69. #define    inbuf    uimg.src
  70. #define    obuf    uimg.dest
  71. #define    frm    uimg.frames
  72. #define    cln    uimg.width
  73. #define    row    uimg.height
  74.  
  75. #ifdef    _DEBUG_
  76. extern    int    debug;
  77. #endif
  78. bool    verbose, spectrum=True, fflip, flip, selfscale, v2dtoc;
  79.  
  80.  
  81. main(argc, argv)
  82. int    argc;
  83. char*    argv[];
  84. {
  85. int    i, f, fsize, dimen1len, vfsize;
  86. double    *dbuf, *vbuf, *cvt, maxi=0, scale=0, divident=1.;
  87.  
  88. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  89. uimg.pxl_out = 1;
  90.  
  91. for (i=1; i<argc; i++)
  92.     if (strcmp(argv[i], "-v")==0)        verbose++;
  93.     else    if (!strncmp(argv[i], "-m", 2) && i+1 < argc)
  94.         maxi = atof(argv[++i]);
  95.     else    if (!strncmp(argv[i], "-di", 3) && i+1 < argc)
  96.         divident = atof(argv[++i]);
  97.     else    if (!strncmp(argv[i], "-C", 2))
  98.         v2dtoc++;
  99.     else    if (!strncmp(argv[i], "-F", 2))
  100.         uimg.pxl_out = sizeof(float);
  101.     else    if (!strncmp(argv[i], "-ff", 3))
  102.         flip = ++fflip;
  103.     else    if (!strncmp(argv[i], "-fl", 3))
  104.         flip++;
  105.     else    if (!strncmp(argv[i], "-r", 2))
  106.         spectrum=0;
  107.     else    if (!strncmp(argv[i], "-s", 2))
  108.         selfscale++;
  109.     else    if (!strcmp(argv[i], "-n"))
  110.         scale = 1.;
  111.     else    if ((in_fp=freopen(argv[i], "rb", stdin))==0)
  112. info:        usage_n_options(usage, i, argv[i]);
  113.  
  114. io_test(fileno(in_fp), goto    info);
  115.  
  116. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  117.  
  118. if (uimg.in_form<IFMT_COMPLEX || uimg.in_form>IFMT_DBLCOM  &&
  119.     (uimg.in_form < IFMT_VFFT3D || uimg.in_form > IFMT_DVVFFT3D))
  120.     syserr("not in double, complex or vfft format");
  121.  
  122. if (uimg.pxl_out == 1)
  123.     uimg.o_form = IFMT_BYTE;
  124. else    uimg.o_form = IFMT_FLOAT;
  125. if (v2dtoc)
  126.     if (uimg.in_form < IFMT_VFFT3D && uimg.in_form > IFMT_DVFFT2D)
  127.     v2dtoc = False;
  128.     else
  129.     mesg("VFFT to Complex FFT\n"),
  130.     uimg.o_form = IFMT_COMPLEX,
  131.     uimg.pxl_out = 8<<(uimg.in_form>=IFMT_DVFFT3D);
  132.  
  133. dimen1len = (cln>>1) + 1;
  134. fsize = row * cln;
  135. if (uimg.in_form > IFMT_DBLCOM)
  136.     vfsize = row * dimen1len;
  137. else    vfsize = fsize;
  138.  
  139. dbuf = nzalloc(vfsize, uimg.pxl_in, "dbuf");
  140. obuf = nzalloc(fsize, uimg.pxl_out, "obuf");
  141. if (vfsize != fsize || flip && uimg.in_form < IFMT_VFFT3D)
  142.     vbuf = nzalloc(fsize, sizeof(*vbuf), "vbuf");
  143.  
  144. (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  145.  
  146. if (v2dtoc)    /* 2D VFFT to complex    */
  147.     for (f=0; f<frm; f++)    {
  148.     int    r_len = row+1>>1, c_len = cln+1>>1;
  149.  
  150.     i=upread(dbuf, uimg.pxl_in, vfsize, in_fp);
  151.     if (i != vfsize)
  152.         syserr("[%d] read %d", vfsize, i);
  153.     if (uimg.in_form==IFMT_VFFT2D || uimg.in_form==IFMT_VFFT3D){
  154.     register float    vscale=scale;
  155.     register COMPLEX    *c_in = PrtCAST dbuf, *cp=obuf;
  156.  
  157.         if (!vscale)    vscale /= fsize * (uimg.in_form & 1 ? frm : 1);
  158.  
  159.         cp->re = c_in->re * vscale;
  160.         cp->im = c_in->im * vscale;
  161.         for (i=1; i<c_len; i++)    {
  162.         cp[cln-i].re = cp[i].re = c_in[i].re * vscale;
  163.         cp[cln-i].im = -(cp[i].im = c_in[i].im * vscale);
  164.         }
  165.         if ((cln & 1)==0){
  166.         cp[i].re = c_in[i].re * vscale;
  167.         cp[i].im = c_in[i].im * vscale;
  168.         }
  169.  
  170.         for (i=1; i<r_len; i++) {    /* dis for conjugate symmetry    */
  171.         register int    j=0, dis=(row-(i<<1))*dimen1len,
  172.             dis1=cln*(row-(i<<1)), dis0=dis1+cln;
  173.         cp += cln;
  174.         c_in += dimen1len;
  175.         cp->re = c_in[j].re * vscale;
  176.         cp->im = c_in[j].im * vscale;
  177.         cp[dis1].re = c_in[dis].re * vscale;
  178.         cp[dis1].im = c_in[dis].im * vscale;
  179.         for (j=1; j<c_len;j++)    {
  180.             cp[dis0-j].re = cp[j].re = c_in[j].re * vscale;
  181.             cp[dis0-j].im = -(cp[j].im = c_in[j].im * vscale);
  182.             cp[cln-j].re = cp[dis1+j].re = c_in[dis+j].re * vscale;
  183.             cp[cln-j].im = -(cp[dis1+j].im = c_in[dis+j].im * vscale);
  184.         }
  185.         if ((cln & 1)==0) {
  186.             cp[j].re = c_in[j].re * vscale;
  187.             cp[j].im = c_in[j].im * vscale;
  188.             cp[dis1+j].re = c_in[dis+j].re * vscale;
  189.             cp[dis1+j].im = c_in[dis+j].im * vscale;
  190.         }
  191.         }
  192.         if ((row & 1)==0)    {
  193.         cp += cln;
  194.         c_in += dimen1len;
  195.         cp->re = c_in->re * vscale;
  196.         cp->im = c_in->im * vscale;
  197.         for (i=1; i<c_len; i++)    {
  198.             cp[cln-i].re = cp[i].re = c_in[i].re * vscale;
  199.             cp[cln-i].im = -(cp[i].im = c_in[i].im * vscale);
  200.         }
  201.         if ((cln & 1)==0){
  202.             cp[i].re = c_in[i].re * vscale;
  203.             cp[i].im = c_in[i].im * vscale;
  204.         }
  205.         }
  206.     }
  207.     else    {
  208.     register double    *ibp = dbuf, vscale=scale;
  209.  
  210.         if (!vscale)    vscale /= fsize * (uimg.in_form & 1 ? frm : 1);
  211.         for (i=0; i<row; i++){
  212.         register int    j;
  213.         register DBCOMPLEX    *cp=obuf;
  214.         cp += i*cln;
  215.         cp->re = *ibp++ * vscale;
  216.         cp->im = *ibp++ * vscale;
  217.         for (j=1; j<cln+1>>1;j++){
  218.             cp[cln-j].re = cp[j].re = *ibp++ * vscale;
  219.             cp[cln-j].im = -(cp[j].im = *ibp++ * vscale);
  220.         }
  221.         if ((cln & 1)==0){
  222.             cp[cln>>1].re = *ibp++ * vscale;
  223.             cp[cln>>1].im = *ibp++ * vscale;
  224.         }
  225.         }
  226.     }
  227.     if ((i=fwrite(obuf, uimg.pxl_out, fsize, out_fp)) != fsize)
  228.         syserr("[%d] write %d", fsize, i);
  229. }
  230. else    for (f=0; f<frm; f++){
  231.     register double    *dp;
  232.  
  233.     cvt = dbuf;
  234.     if ((i=upread(dbuf, uimg.pxl_in, vfsize, in_fp)) != vfsize)
  235.         syserr("[%d] Read %d", vfsize, i);
  236.     if (uimg.in_form==IFMT_COMPLEX)
  237.         complex_to_d(dbuf, vfsize);
  238.     else if (uimg.in_form == IFMT_DBLCOM)
  239.         dcomplex_to_d(dbuf, vfsize);
  240.     else if (uimg.in_form > IFMT_DBLCOM)
  241.         vfft_handle(uimg.in_form, dbuf, cvt=vbuf, row, cln, dimen1len);
  242.     if (uimg.in_form < IFMT_VFFT3D)
  243.         if (fflip)
  244.         Dcom_fflip(dbuf, cvt=vbuf, row, cln);
  245.         else if (flip)
  246.         Dcom_flip(dbuf, cvt=vbuf, row, cln);
  247.  
  248.     dp = cvt;
  249.  
  250.     if (uimg.o_form == IFMT_BYTE){
  251.     register byte    *bp = obuf;
  252.     double    imin;
  253.  
  254.         if (scale != 1. && f==0 || selfscale || scale > 1.7e308){
  255.         register double    min=1.7e308, max=1.7e-308, value;
  256.         for (i=0; i<fsize; i++){
  257.             value = dp[i];
  258.             if (value < min)    min = value;
  259.             else if (value > max)    max = value;
  260.         }
  261.         if (verbose || uimg.in_form > IFMT_DBLCOM)
  262.             message("%s: min_in=%f, max_in=%f\n", Progname, min, max);
  263.         if (maxi)
  264.             max = maxi;
  265.         else    max /= divident;
  266.         scale = 255. / (max - min);
  267.         imin = min;
  268.         if (verbose || uimg.in_form > IFMT_DBLCOM)
  269.             message("%s: max_in set to=%4.4f, scale=%f\n",
  270.             Progname, max, scale);
  271.         }
  272.         {
  273.         register double    bscale = scale, min = imin;
  274.         if (bscale > 1.7e308)
  275.             bscale = 1.;
  276.         for (i=0; i<fsize; i++){
  277.         register int tmp = bscale * (dp[i] - min);
  278.             if (tmp > 255)
  279.                 bp[i] = 255;
  280.             else    bp[i] = tmp;
  281.         }
  282.         }
  283.     }
  284.     else{
  285.     register float    *fp = obuf;
  286.         if (uimg.in_form > IFMT_DBLCOM){
  287.         register double    fscale = 1. / fsize;
  288.         for (i=0; i<fsize; i++)
  289.             *fp++ = fscale * *dp++;
  290.         }
  291.         else for (i=0; i<fsize; i++)
  292.         *fp++ = *dp++;
  293.     }
  294.     if ((i=fwrite(obuf, uimg.pxl_out, fsize, out_fp)) != fsize)
  295.         syserr("[%d] write %d", fsize, i);
  296. }
  297. }
  298.  
  299.  
  300. /*    float complex to double    */
  301. complex_to_d(buf, size)
  302. VType    *buf;
  303. {
  304. register float    *fp = buf;
  305. register double    *op = buf;
  306. register int    i;
  307.  
  308. if (verbose)    mesg("converting complex to double\n");
  309. for (i=0; i<size; i++){
  310. register double    value = *fp * *fp;
  311.     fp++;
  312.     value += *fp * *fp;
  313.     fp++;
  314.     *op++ = sqrt(value);
  315. }
  316. }
  317.  
  318. /*    double complex to real double    */
  319. dcomplex_to_d(buf, size)
  320. register double    *buf;
  321. {
  322. register double    *ibp = buf;
  323. register int    i;
  324.  
  325. if (verbose)    mesg("converting double complex to double\n");
  326. for (i=0; i<size; i++){
  327. register double    value = *ibp * *ibp;
  328.     ibp++;
  329.     value += *ibp * *ibp;
  330.     ibp++;
  331.     *buf++ = sqrt(value);
  332. }
  333. }
  334.  
  335.  
  336. /*===============================================
  337. %    This routine convert vfft to spectrum.    %
  338. %        It is exactly transform.    %
  339. ===============================================*/
  340. vfft_handle(vform, ibuf, o_buf, rows, cols, dimen1len)
  341. double    *ibuf, *o_buf;
  342. {
  343. int    vsize = rows * dimen1len;
  344.  
  345. switch(vform){
  346. case IFMT_VFFT2D:
  347. case IFMT_VFFT3D:    complex_to_d(ibuf, vsize);
  348.     break;
  349. case IFMT_DVFFT2D:
  350. case IFMT_DVFFT3D:    dcomplex_to_d(ibuf, vsize);
  351.     break;
  352. }
  353.  
  354. if (spectrum)    {    /*    VFFT spectrum    */
  355. register int    i, j, dc=dimen1len-(dimen1len&1);
  356. register double    *src=ibuf + (vsize>>1),    *dst=o_buf + dimen1len-1,
  357.     *dstcs=o_buf + rows*cols - (dimen1len + (cols&1));
  358.     /* like in v2dtoc, we can use integer dis instead of using *dstcs
  359.         for computing conjugate symmetry distance */
  360.     for (i=0; i<rows>>1; i++)    {
  361.     *dst = *dstcs = *src++;
  362.     for (j=1; j<dc; j++)
  363.         dstcs[-j] = dst[j] = *src++;
  364.     if (dimen1len & 1)    /* throw away the Nyquist since just display */
  365.         src++;
  366.     dst += cols;    dstcs -= cols;
  367.     }
  368.  
  369.     for (src=ibuf; i<rows; i++)    {
  370.     *dst = *dstcs = *src++;
  371.     for (j=1; j<dc; j++)
  372.         dstcs[-j] = dst[j] = *src++;
  373.     if (dimen1len & 1)
  374.         src++;
  375.     dst += cols;    dstcs -= cols;
  376.     }
  377. }
  378. else    {
  379. register int    i, j, dis=cols*rows;    /* conjugately symmetric distance */
  380. register double    *src=ibuf,    *dst=o_buf;
  381.     for (i=0; i<rows; i++)    {
  382.     *dst = *src++;
  383.     for (j=1; j<dimen1len; j++)
  384.         dst[dis-j] = dst[j] = *src++;
  385.     dst += cols;    dis -= cols << 1;
  386.     }
  387. }
  388. }
  389.  
  390. /*    Double Complex Flip Quadrant    */
  391. Dcom_flip(ibuf, o_buf, rows, cols)
  392. double    *ibuf, *o_buf;
  393. {
  394. int    hcols=cols>>1;
  395. register int    i, j;
  396. register double    *src=ibuf + (rows*cols>>1),    *dst=o_buf + hcols;
  397.     for (i=0; i<rows>>1; i++){
  398.     for (j=0; j<hcols; j++)
  399.         dst[-j-1] = dst[j] = *src++;
  400.     dst += cols;    src += hcols;
  401.     }
  402.     for (src=ibuf; i<rows; i++){
  403.     for (j=0; j<hcols; j++)
  404.         dst[-j-1] = dst[j] = *src++;
  405.     dst += cols;    src += hcols;
  406.     }
  407. }
  408.  
  409. /*===============================================
  410. %    left half is one line higher than    %
  411. %    right half. It does not effect view,    %
  412. %    but may not use for real application    %
  413. ===============================================*/
  414. Dcom_fflip(ibuf, o_buf, rows, cols)
  415. double    *ibuf, *o_buf;
  416. {
  417. register int    i;
  418. register double    *src=ibuf + (cols*(rows+1) >> 1), *dst = o_buf;
  419.  
  420. for (i=cols*(rows-1) >> 1; i--;)
  421.     *dst++ = *src++;
  422. src = ibuf;
  423. for (i=cols*(rows+1) >> 1; i--;)
  424.     *dst++ = *src++;
  425. }
  426.