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