home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11.lha / ccs-lib / lib / fitsread.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-14  |  8.7 KB  |  302 lines
  1. /*
  2. %    FITSREAD . C
  3. @
  4. @    read FITS images
  5. @
  6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  7.  
  8. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  9. Permission is granted to reproduce this software for non-commercial
  10. purposes provided that this notice is left intact.
  11.  
  12. It is acknowledged that the U.S. Government has rights to this software
  13. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  14. and the University of California.
  15.  
  16. This software is provided as a professional and academic contribution
  17. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  18. with no warranties of any kind whatsoever, no support, no promise of
  19. updates, or printed documentation. By using this software, you
  20. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  21. University of California shall have no liability with respect to the
  22. infringement of other copyrights by any part of this software.
  23.  
  24. For further information about this notice, contact William Johnston,
  25. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  26. (wejohnston@lbl.gov)
  27.  
  28. For further information about this software, contact:
  29.     Jin Guojun
  30.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  31.     g_jin@lbl.gov
  32.  
  33. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  34. @
  35. @ AUTHOR:    Guojun Jin - LBL    9/1/91
  36. */
  37.  
  38. #include "header.def"
  39. #include "fits2hip.h"
  40. #include <math.h>
  41.  
  42. #ifndef    BufSize
  43. #define    BufSize    4096
  44. #endif    114030 = 1BD6E for unix
  45.  
  46. /************************************************************************
  47. *    uncompress is based on the previously difference algorithm.    *
  48. *    If input is not equal to -128, which means the difference    *
  49. *    is less than 255 (maximum byte value).    So, we just simple    *
  50. *    add current  value with pre byte value    to get current value.    *
  51. *    Otherwise, goto Auxiliary buffer to pick the long value.    *
  52. *    J is Auxiliary buffer pointer. I is input buffer pointer.    *
  53. *    After decompression, if J is not equal to Naux, then there    *
  54. *    must be some thing wrong. So, please check J value.        *
  55. *    Warning:                            *
  56. *        This uncompress will not properly work on IBM/PC for    *
  57. *    large file. Since it must store all date blocks in main memory.    *
  58. ************************************************************************/
  59.  
  60. char    fits_uncompress(img, FTy)
  61. U_IMAGE    *img;
  62. bool    FTy;
  63. {
  64. VType    *unc_ibuf;
  65. float    *unc_cnvt=NULL;
  66. long    prev, next, j=0, length, fsize=img->height*img->width,
  67.     NP=fsize*img->frames;
  68. register FITSType*    aux = zalloc(Naux, 2);
  69. register char*    p;
  70. register float*    fp;
  71. register MType    i;
  72. #ifdef    _DEBUG_
  73. FITSType    maxval=MaxShort, minval = MinShort;
  74. #endif
  75.  
  76. p = unc_ibuf = nzalloc(img->frames, fsize, "unc_in");
  77. if (img->pxl_out != sizeof(*unc_cnvt))
  78.     unc_cnvt = nzalloc(NP, sizeof(*unc_cnvt), "unc_cnvt");
  79. else    unc_cnvt = img->src;
  80.  
  81.     if (FORTRAN){
  82.         for (length=0; length<NP;)    /* reading difference data */
  83.         length += read_var(p+length, img->IN_FP, hostype, FTy);
  84.  
  85.         for (length=0; length < Naux<<1L;)    /* reading auxiliary data */
  86.         length += read_var(aux + length, img->IN_FP, hostype, FTy);
  87.     }
  88.     else{
  89.         message("attention: Suppose not to be a FORTRAN File\n");
  90.         fread(p, img->frames, fsize, img->IN_FP);
  91.         fread(aux, Naux, 2, img->IN_FP);
  92.     }
  93.  
  94. if (((hostype==2 || hostype==5) && FTy != 'u') || (hostype != 5 && FTy == 'u'))
  95.     for (i=0; i < Naux; i++)
  96.     *(aux + i) = ShortSwap(*(aux + i));
  97. else    msg("No swaping, FTy = %c\n", FTy);
  98.  
  99.     prev = *aux;
  100.     DEBUGMESSAGE("First Aux is %d,    first value = %d\n", prev, prev + *p);
  101.  
  102.     for (i=NP, fp=unc_cnvt; i--; p++){
  103.         if (*p != -128)
  104.         next = prev + *p;
  105.         else{
  106.         if    (j >= Naux)
  107.             return    NULL;
  108.         next = aux[++j];
  109.         }
  110.         *fp++ = prev = next;
  111.  
  112. #ifdef    _DEBUG_
  113.         if (maxval < prev)    maxval = prev;
  114.         if (minval > prev)    minval = prev;
  115. #endif
  116.     }
  117.     free(aux);
  118.     free(unc_ibuf);
  119. #ifdef    _DEBUG_
  120.     msg("DINF: Min_value = %4ld,    Max_value = %ld\n", minval, maxval);
  121. #endif
  122.  
  123.     if (Scale && Scale != 1)
  124.         for (i=0, fp=unc_cnvt; i < NP; i++, fp++)
  125.         *fp = *fp * Scale + Bzero;
  126.     fp = unc_cnvt;
  127.     switch(img->o_form){
  128.     case IFMT_BYTE:
  129.         for (i=0, p=img->src; i < NP; i++, p++, fp++){
  130.             *fp -= ZLevel;
  131.             if (*fp < 256)
  132.                 *(byte*)p = *fp;
  133.             else    *p = 255;
  134.         }
  135.         break;
  136.     case IFMT_SHORT:
  137.         for (i=0, aux=img->src; i < NP; i++)
  138.             *aux++ = *fp++ - ZLevel;
  139.         break;
  140.     case IFMT_LONG:
  141.         for (i=0; i < NP; i++)
  142.         ((long*)img->src)[i] = *fp++ - ZLevel;
  143.     }
  144.     if (img->pxl_out != sizeof(*unc_cnvt))
  145.         free(unc_cnvt);
  146.     if (j != Naux - 1)
  147.         msg("diff Naux(%d) = %d, F %d\n", Naux, j, feof(img->IN_FP));
  148. return    (Naux - 1 - j);
  149. }
  150.  
  151.  
  152. /************************************************************************
  153. *    This procedure only swap and convert data from fits format to    *
  154. *    regular data format. You may need to use mainpeak to get    *
  155. *    proper picture.    (called by transf_data)                *
  156. ************************************************************************/
  157.  
  158. float    fits_convertdata(img, ibuf, length, obuf)
  159. U_IMAGE    *img;
  160. long    *ibuf, length;
  161. register VType    *obuf;
  162. {
  163. float    *buf=nzalloc((MType)sizeof(float), (MType)BufSize, "fits_cnvt");
  164. register FITSType*    sp;
  165. register float    *fp = buf, mean=0;
  166. register long    i;
  167. register float    coef = Bzero - Seeing * ZLevel, Fscale = Bscale;
  168.  
  169. if (((hostype==2 || hostype==5) && FTy != 'u') ||    /* u & v => t    */
  170.     (hostype != 2 && hostype != 5 && FTy == 'u'))
  171.     switch(img->pxl_in)    {
  172.     case 2:
  173.         for (i=0, sp=PrtCAST ibuf; i < length; i++)
  174.         {    *sp = ShortSwap(*sp);
  175.             fp[i] = *sp++ * Fscale + coef;
  176.         }
  177.     break;
  178.     case 4:
  179.         for (i=0; i < length; i++)
  180.         {    *(ibuf + i) = LongSwap(*(ibuf + i));
  181.             fp[i] = *(ibuf + i) * Fscale + coef;
  182.         }
  183.     }
  184. else    switch(img->pxl_in) {
  185.     case 2:    for (i=0, sp=PrtCAST ibuf; i < length; i++)
  186.             fp[i] = *sp++ * Fscale + coef;
  187.         break;
  188.     case 4:    for (i=0; i < length; i++)
  189.             fp[i] = *(ibuf + i) * Fscale + coef;
  190.     }
  191.  
  192. switch (img->o_form)    {
  193. case IFMT_BYTE:
  194.     for (i=0; i < length; mean += fp[i++])
  195.         if (fp[i] > 255)
  196.             *((char*)obuf + i) = 255;
  197.         else if (fp[i] < 0)
  198.             *((char*)obuf + i) = 0;
  199.         else    *((char*)obuf + i) = fp[i];
  200.     break;
  201. case IFMT_LONG:
  202.     for (i=0; i < length; mean += fp[i++])
  203.         ((long*)obuf)[i] = fp[i];
  204.     break;
  205. case IFMT_FLOAT:
  206.         for (i=0; i < length; mean += fp[i++]);
  207.         memcpy(obuf, fp, length*sizeof(*fp));
  208.     break;
  209. case IFMT_DOUBLE:
  210.     for (i=0; i < length; mean += *fp++, i++)
  211.         ((double*)obuf)[i] = *fp;
  212.     break;
  213. case IFMT_SHORT:
  214. default:
  215.     for (i=0, sp=obuf; i < length; mean += *fp++, i++, sp++)
  216.         if (*fp > MaxShort)
  217.             *sp = MaxShort;
  218.         else if (*fp < MinShort)
  219.             *sp = MinShort;
  220.         else    *sp = *fp;
  221. }
  222. free(buf);
  223. return    mean;
  224. }
  225.  
  226. /************************************************************************
  227. *    There is no decompression. It only transform fits image data    *
  228. *    to hips image data. It can be used on variant machines.        *
  229. ************************************************************************/
  230. fits_transf_data(img, FTy, wflag)
  231. U_IMAGE    *img;
  232. bool    FTy, wflag;
  233. {
  234. MType    block, length, pixels,
  235.     isize=img->pxl_in*img->frames*img->height*img->width;
  236. float    sumean=0;
  237. VType    *t_ibuf = nzalloc(img->pxl_in, BufSize, "transf_in"),
  238.     *t_cnvt = nzalloc(img->pxl_out, BufSize, "transf_cnvt");
  239.  
  240. #ifdef    _DEBUG_
  241. msg("SEEing = %f, CrVal1 = %f\n", Seeing, Crval1);
  242. #endif
  243.  
  244. if (FORTRAN)
  245.     for (length=block=0; length < isize;)    {
  246.     pixels = read_var(t_ibuf, img->IN_FP, hostype, FTy);
  247.     if (pixels < 0 || feof(img->IN_FP)) {
  248.         mesg("Warning!    may have errors\n");
  249.         break;
  250.     }
  251.     length += pixels;
  252.     pixels >>= ByteShift;        /* get pixel number    */
  253.     sumean += fits_convertdata(img, t_ibuf, pixels, t_cnvt);
  254.     if (wflag)    fwrite(t_cnvt, oPBSize, pixels, img->OUT_FP);
  255.     else    memcpy((char*)img->src+length, t_cnvt, pixels<<ByteShift);
  256.     }
  257. else {    mesg("Attention: Suppose not to be a FORTRAN File\n");
  258.  
  259.     for (length=0, pixels=BufSize; length < isize && pixels==BufSize;
  260.         length += pixels<<ByteShift)
  261.     {
  262.         pixels = fread(t_ibuf, 1<<ByteShift, BufSize, img->IN_FP);
  263.         sumean += fits_convertdata(img, t_ibuf, pixels, t_cnvt);
  264.         if (wflag)
  265.             fwrite(t_cnvt, oPBSize, pixels, img->OUT_FP);
  266.         else    memcpy((char*)img->src+img->pxl_out*(length>>ByteShift),
  267.                 t_cnvt, pixels*img->pxl_out);
  268.     }
  269. }
  270. msg("Mean Value = %lf\n", sumean/isize);
  271. free(t_ibuf);
  272. free(t_cnvt);
  273. return    (int) length - (block << (FTy != 'v')) * FFCL - isize;
  274. }
  275.  
  276.  
  277. read_fits_image(img, Fty, dcmp, wflag)
  278. U_IMAGE    *img;
  279. bool    Fty, dcmp, wflag;    /* write data to a file */
  280. {
  281. int    errcode;
  282.  
  283. if (dcmp)    {    /* decompress */
  284.     verify_buffer_size(&img->src, img->frames * img->pxl_out,
  285.             img->height * img->width, "unc_src");
  286.     if (errcode=fits_uncompress(img, Fty))
  287.         msg("%d error(s) in decompressing\n", errcode);
  288.     else    mesg("end of uncompressing\n");
  289.     if (wflag)
  290.         update(img->src, img->pxl_out, (MType)img->height * img->width,
  291.             img->OUT_FP);
  292. }
  293. else    {
  294.     if (!wflag)
  295.         verify_buffer_size(&img->src, img->frames * img->pxl_out,
  296.             img->height * img->width, "transf_src");
  297.     errcode = fits_transf_data(img, Fty, wflag);
  298.     msg("error(%d)    end of FITS transfer\n", errcode);
  299. }
  300. return    errcode;
  301. }
  302.