home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / sun / lib2 / dbvfft3d.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-04  |  6.0 KB  |  260 lines
  1. /*
  2. %    DBVFFT_3D . 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. %    The 1st dimension is shrinked to (m/2 + 1), buf
  34. %    the 2nd & 3rd dimensions are doubled by complex format.
  35. %
  36. % AUTHOR:    Jin Guojun - LBL    5/1/91
  37. */
  38.  
  39. #include    "complex.h"
  40. #include    "imagedef.h"
  41.  
  42. DBCOMPLEX    *DBcmpx_in, *DBcmpx_out;
  43.  
  44. void    DBwork_space_init(max_len)
  45. {
  46. static MType    space_len=0;
  47.  
  48. if (space_len < max_len){
  49.     space_len = max_len;
  50.     if (DBcmpx_in)    free(DBcmpx_in);
  51.     if (DBcmpx_out)    free(DBcmpx_out);
  52.     DBcmpx_in = (DBCOMPLEX*) nzalloc(space_len, (MType)sizeof(*DBcmpx_in), "cmx_in");
  53.     DBcmpx_out = (DBCOMPLEX*)nzalloc(space_len, (MType)sizeof(*DBcmpx_out), "cmx_out");
  54. }
  55. }
  56.  
  57. double    *DBvfft2d(src, x, y, buf)
  58. VType    *src;
  59. unsigned x, y;
  60. VType    *buf;
  61. {
  62. double    *in = (double *) src, *out = (double *) buf;
  63. register unsigned    i, loop, n;
  64. register DBCOMPLEX    *cmptr;
  65.  
  66. if (!x || load_DBw(x) == -1)    return    in;
  67.  
  68. if (x>1){
  69.     for (n=x, loop=0; loop<y; loop++){
  70.     register double*    dp = in + loop * n;
  71.     cmptr = DBcmpx_in;
  72.     for (i = 0; i < n; i++) {
  73.         c_re (cmptr [i]) = dp[i];
  74.         c_im (cmptr [i]) = 0;
  75.     }
  76.  
  77.     DBFourier(cmptr, n, DBcmpx_out);
  78.  
  79.     {
  80.     register int    len = ((n>>1) + 1) << 1;
  81.     dp = out + loop * len; /* load x/2 + 1  COMPLEXs */
  82.     memcpy(dp, DBcmpx_out, len * sizeof(*dp));
  83.     }
  84. #ifdef    _DEBUG_
  85.     printcmpl(cmptr, n, loop, "FX");
  86. #endif
  87.     }
  88.     in = (double*)buf;
  89. }
  90.  
  91. if (y<2 || load_DBw(y) == -1)    return    in;
  92. #ifdef    _DEBUG_
  93. printsamp(in, x*y);    /* here, in=out=buf */
  94. #endif
  95.  
  96. {
  97. register int    dimen1len = (x>>1) + 1;/* 1st dimension is shorted to x/2 + 1 */
  98.  
  99.     for (n=y, loop=dimen1len; loop--;){
  100.     register DBCOMPLEX *cp = (DBCOMPLEX*)in + loop; /* carefully convert */
  101.         cmptr = DBcmpx_in;
  102.         for (i = 0; i < n; i++, cp+=dimen1len)
  103.             cmptr[i] = *cp;
  104.  
  105.         DBFourier(cmptr, n, DBcmpx_out);
  106.  
  107.         cp = (DBCOMPLEX*)out + loop;
  108.         cmptr = DBcmpx_out;
  109.         for (i=0; i < n; i++, cp+=dimen1len)
  110.             *cp = cmptr[i];
  111.     }
  112. }
  113. return    out;
  114. }
  115.  
  116.  
  117. double    *DBvrft2d(src, x, y, buf)
  118. VType    *src;
  119. unsigned x, y;
  120. VType    *buf;
  121. {
  122. double    *in=(double *) src, *out=(double *) buf;
  123. register unsigned    i, loop, n;
  124. register DBCOMPLEX    *cmptr;
  125. register int    dimension1 = (x>>1) + 1; /* for input only */
  126.  
  127. if (!x || load_DBw(y) == -1)    return    in;
  128.  
  129. if (y>1) {
  130.  
  131.    for (n=y, loop=dimension1; loop--;){
  132.    register DBCOMPLEX    *cp = (DBCOMPLEX*)in + loop;
  133.     cmptr = DBcmpx_in;
  134.     for (i = 0; i < n; i++, cp+=dimension1)
  135.         cmptr[i] = *cp;
  136.  
  137.     DBFourier(cmptr, n, DBcmpx_out);
  138.  
  139.     cp = (DBCOMPLEX*)in + loop;    /* back to src */
  140.     cmptr = DBcmpx_out;
  141.     for (i=0; i < n; i++, cp+=dimension1)
  142.         *cp = cmptr[i];
  143. #ifdef    _DEBUG_
  144.     printcmpl(cmptr, n, loop, "RY");
  145. #endif
  146.     }
  147. }
  148.  
  149. if (x<2 || load_DBw(x) == -1)    return    in;
  150. #ifdef    _DEBUG_
  151. printsamp(in, x*y);    /* here, in=src, out=buf */
  152. #endif
  153.  
  154. for (n=x, loop=0; loop<y; loop++){
  155. register DBCOMPLEX*    cp = (DBCOMPLEX*)in + loop*dimension1;
  156.     cmptr = DBcmpx_in;
  157.     cmptr [0] = cp[0];        /* dc */
  158.     for (i=(n+1 >> 1); --i;) {    /* conj. symm. harmonischen */
  159.         cmptr[i] = cp[i];
  160.         c_re(cmptr[n-i]) = c_re(cp[i]);
  161.         c_im(cmptr[n-i]) = -c_im(cp[i]);
  162.     }
  163.     if ((n & 1) == 0)        /* Nyquist */
  164.         cmptr[n >> 1] = cp[n >> 1];
  165.  
  166.     DBFourier(cmptr, n, DBcmpx_out);
  167.  
  168.     {
  169.     register double    *dp = out + loop*n;
  170.         for (i=0, cmptr=DBcmpx_out; i < n; i++)
  171.         dp[i] = c_re (cmptr [i]);
  172.     }
  173. #ifdef    _DEBUG_
  174.     printcmpl(cmptr, n, loop, "RX");
  175. #endif
  176. }
  177. #ifdef    _DEBUG_
  178. printsamp(out, x*y);
  179. #endif
  180. return    out;
  181. }
  182.  
  183. VType    *DBvfft3d(src, x, y, z, dst)
  184. double    *src;
  185. DBCOMPLEX    *dst;
  186. {
  187. register int    i, dimen1len = (x>>1) + 1, dimen2size=dimen1len * y;
  188. int    fsize=x*y;
  189.  
  190. DBwork_space_init(MAX(MAX(x, y), z));
  191. for (i=0; i<z; i++)
  192.     DBvfft2d(src+i*fsize, x, y, dst+i*dimen2size);
  193. if (z>1){
  194. int    r;
  195. register DBCOMPLEX    *cmptr, *cp;
  196. register int    c, n=z;
  197.  
  198. load_DBw(z);
  199.  
  200. for (r=0; r<y; r++)
  201.     for (c=0; c<dimen1len; c++){
  202.     cp = dst + r*dimen1len+c;
  203.     cmptr = DBcmpx_in;
  204.     for (i = 0; i < n; i++, cp+=dimen2size)
  205.         cmptr[i] = *cp;
  206.  
  207.     DBFourier(cmptr, n, DBcmpx_out);
  208.  
  209.     cp = dst + r*dimen1len+c;    /* store back to dst */
  210.     cmptr = DBcmpx_out;
  211.     for (i = 0; i < n; i++, cp+=dimen2size)
  212.         *cp = cmptr[i];
  213. #ifdef    _DEBUG_
  214.     printcmpl(cmptr, n, r, "FZ");
  215. #endif
  216.     }
  217. }
  218. return    (VType*)dst;
  219. }
  220.  
  221. VType    *DBvrft3d(src, x, y, z, dst)
  222. DBCOMPLEX    *src;
  223. double    *dst;
  224. {
  225. register int    i, dimen1len=(x>>1)+1, dimen2size=dimen1len*y;
  226. int    fsize=x*y;
  227.  
  228. if (z>1){
  229. int    r;
  230. register DBCOMPLEX    *cmptr, *cp;
  231. register int    c, n=z;
  232.  
  233. DBwork_space_init(MAX(MAX(x, y), z));
  234. load_DBw(n);
  235.  
  236. for (r=0; r<y; r++)
  237.     for (c=0; c<dimen1len; c++){
  238.     cp = src + r*dimen1len+c;
  239.     cmptr = DBcmpx_in;
  240.     for (i = 0; i < n; i++, cp+=dimen2size) {
  241.         c_re (cmptr [i]) = c_re(*cp);
  242.         c_im (cmptr [i]) = -c_im(*cp);
  243.     }
  244.  
  245.     DBFourier(cmptr, n, DBcmpx_out);
  246.  
  247.     cp = src + r*dimen1len+c;    /* store back to inbuf */
  248.     cmptr = DBcmpx_out;
  249.     for (i = 0; i < n; i++, cp+=dimen2size)
  250.         *cp = cmptr[i];
  251. #ifdef    _DEBUG_
  252.     printcmpl(cmptr, n, r, "FZ");
  253. #endif
  254.     }
  255. }
  256. for (i=0; i<z; i++)
  257.     vrft2d(src + i*dimen2size, x, y, dst + i*fsize);
  258. return    (VType*)dst;
  259. }
  260.