git.oblomov.eu Git - zanmar/loader_fits/blob - loader_fits.c

   1 #include <loader_common.h>
   2 #include <fitsio.h>
   3
   4 //TODO: add a color table.
   5 //invert gray scale, Check if feh or qiv can do this.
   6 // add option to check for datamin and datamax on the header
   7 // make it more robust. datamax > datamin
   8 // write the fits writer!
   9
  10 #define SWAP32(x) (x) = \
  11 ((((x) & 0x000000ff ) << 24) |\
  12  (((x) & 0x0000ff00 ) << 8) |\
  13  (((x) & 0x00ff0000 ) >> 8) |\
  14  (((x) & 0xff000000 ) >> 24))
  15
  16 char
  17 load(ImlibImage * im, ImlibProgressFunction progress,
  18      char progress_granularity, char immediate_load)
  19 {
  20     fprintf(stderr, "fits load function\n");
  21
  22    int                 w = 0, h = 0, alpha = 0;
  23    //FILE               *f;
  24    fitsfile            *f;
  25    int                 status = 0;
  26
  27    if (im->data)
  28       return 0;
  29    //f = fopen(im->real_file, "rb");
  30    if ( fits_open_file( &f, im->real_file, READONLY, &status) ) {
  31          return 0;
  32    }
  33
  34    //if (!f)
  35    //   return 0;
  36
  37    /* header */
  38    {
  39       long                naxes[2];
  40       int                 nfound;
  41       //char                buf[256], buf2[256];
  42
  43       /* read the NAXIS1 and NAXIS2 keyword to get image size */
  44       if ( fits_read_keys_lng( f, "NAXIS", 1, 2, naxes, &nfound, &status ) ){
  45            fits_close_file(f, &status);
  46            return 0;
  47       }
  48       w = naxes[ 0 ];
  49       h = naxes[ 1 ];
  50
  51       /*
  52       buf[0] = '\0';
  53       if (!fgets(buf, 255, f))
  54         {
  55            fclose(f);
  56            return 0;
  57         }
  58       buf2[0] = '\0';
  59       sscanf(buf, "%s %i %i %i", buf2, &w, &h, &alpha);
  60       if (strcmp(buf2, "ARGB"))
  61         {
  62            fclose(f);
  63            return 0;
  64         }
  65         */
  66       if (!IMAGE_DIMENSIONS_OK(w, h))
  67         {
  68            fits_close_file(f, &status);
  69            //fclose(f);
  70            return 0;
  71         }
  72       im->w = w;
  73       im->h = h;
  74       if (!im->format)
  75         {
  76            if (alpha)
  77               SET_FLAG(im->flags, F_HAS_ALPHA);
  78            else
  79               UNSET_FLAG(im->flags, F_HAS_ALPHA);
  80            im->format = strdup("fits");
  81         }
  82    }
  83    if (((!im->data) && (im->loader)) || (immediate_load) || (progress))
  84      {
  85         DATA32             *ptr;
  86         DATA8              gray;
  87         int                 y, x, pl = 0;
  88         char                pper = 0;
  89         double             *buffer, datamin, datamax, nullval = 0;
  90         int                anynull, nbuffer;
  91         long               fpixel = 1;
  92         long               *bin, nmax = 0, ndataok, sumpix = 0;
  93         int                nbin, k, kmin, mostfreq = 0;
  94         double             binw;
  95         char               *envdmin, *envdmax;
  96
  97         /* must set the im->data member before callign progress function */
  98         ptr = im->data = malloc(w * h * sizeof(DATA32));
  99         nbuffer = w;
 100         buffer = malloc( nbuffer * sizeof( double ) );
 101
 102         if (!im->data)
 103           {
 104              fits_close_file(f, &status);
 105              //fclose(f);
 106              return 0;
 107           }
 108
 109         envdmin = getenv( "DATAMIN" );
 110         envdmax = getenv( "DATAMAX" );
 111         if( envdmin != NULL && envdmax != NULL ){
 112             datamin = atof( envdmin );
 113             datamax = atof( envdmax );
 114             fprintf(stderr, "ENV:data min %g, max %g\n", datamin, datamax);}
 115         else{
 116
 117             datamin = HUGE_VALF;
 118             datamax = -HUGE_VALF;
 119
 120             //Find min max
 121             for (y = 0; y < h; y++)
 122               {
 123                  fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 124                               buffer, &anynull, &status);
 125                  for( x = 0; x < w; x++ ){
 126                      if( buffer[ x ] > datamax )
 127                          datamax = buffer[ x ];
 128                      if( buffer[ x ] < datamin )
 129                          datamin = buffer[ x ];
 130                  }
 131                  fpixel += nbuffer;
 132               }
 133             fprintf(stderr, "data min %g, max %g\n", datamin, datamax);
 134
 135             nbin = 4 * sqrt( datamax - datamin );
 136             if( nbin > 5000 ) nbin = 1000;
 137             binw = ( datamax - datamin ) / nbin;
 138             fprintf(stderr, "nbin and bin width:%d %g\n", nbin, binw);
 139             bin = calloc( nbin, sizeof( long ) );
 140
 141             //Find histogram
 142             fpixel = 1;
 143             for (y = 0; y < h; y++)
 144               {
 145                  fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 146                               buffer, &anynull, &status);
 147                  for( x = 0; x < w; x++ ){
 148                      for( k = 0; k < nbin; k++ ){
 149                          if( buffer[ x ] >= k * binw &&
 150                              buffer[ x ] < ( k + 1 ) * binw ){
 151                              bin[ k ] += 1;
 152                              break;
 153                          }
 154                      }
 155                  }
 156                  fpixel += nbuffer;
 157               }
 158             //Find most frequent pixel value
 159             for( k = 0; k < nbin; k++ ){
 160                 if( bin[ k ] > nmax ){
 161                     nmax = bin[ k ];
 162                     mostfreq = k;
 163                 }
 164                 fprintf(stderr, "%ld ", bin[ k ] );
 165             }
 166             fprintf(stderr, "mostfreq and nmax = %d %ld ", mostfreq, nmax );
 167
 168             //Set new DATAMIN. look for FWHM on a log scale
 169             nmax = log( nmax ) / 2;
 170             datamin = ( mostfreq + 0.5 ) * binw;
 171             for( k = mostfreq - 1; k >= 0; k-- ){
 172                 if( bin[ k ] > 0 )
 173                     if( log( bin[ k ] ) < nmax ){
 174                         kmin = k + 1;
 175                         datamin = ( (kmin) + 0.5 ) * binw;
 176                         break;
 177                     }
 178             }
 179             if( k == 0 ){   //case were we have one single peak surrounded by 0
 180                kmin = mostfreq;
 181                datamin = ( (kmin) + 0.5 ) * binw;
 182             }
 183
 184             //Set new DATAMAX when we have 99.5% of the pixels accounted for
 185             ndataok = h * w;
 186             for( k = 0; k < kmin; k++)
 187                 ndataok -= bin[ k ];
 188             //ndataok = h * w; //To get values closer to ds9's uncomment this
 189             ndataok = 0.997 * ndataok;
 190             for( k = kmin; k < nbin; k++){
 191                 sumpix += bin[ k ];
 192                 if( sumpix > ndataok )
 193                     break;
 194             }
 195             datamax = ( ( k + 1 ) + 0.5 ) * binw;
 196         }
 197
 198         // use the new datamin datamax values to rescale and create imlib im
 199         ptr += w * h - 1;
 200         fpixel = 1;
 201         for (y = 0; y < h; y++)
 202           {
 203              fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 204                           buffer, &anynull, &status);
 205
 206              for( x = w - 1; x >= 0; x-- ){
 207                  if( buffer[ x ] <= datamin )
 208                      gray = 0;
 209                  else
 210                      if( buffer[ x ] >= datamax )
 211                          gray = 255;
 212                      else
 213                          gray = ceil(( buffer[ x ] - datamin ) / ( datamax - datamin ) * 255 );
 214                  *ptr = (((( (0xff << 8) | gray) << 8) | gray) << 8) | gray ;
 215                  ptr--;
 216              }
 217              fpixel += nbuffer;
 218
 219              if (progress)
 220                {
 221                   char                per;
 222                   int                 l;
 223
 224                   per = (char)((100 * y) / im->h);
 225                   if (((per - pper) >= progress_granularity) ||
 226                       (y == (im->h - 1)))
 227                     {
 228                        l = y - pl;
 229                        if (!progress(im, per, 0, (y - l), im->w, l))
 230                          {
 231                             fits_close_file(f, &status);
 232                             //fclose(f);
 233                             return 2;
 234                          }
 235                        pper = per;
 236                        pl = y;
 237                     }
 238                }
 239           }
 240             fprintf(stderr, "data min %g, max %g\n", datamin, datamax);
 241      }
 242    //fclose(f);
 243    fits_close_file(f, &status);
 244    return 1;
 245 }
 246
 247 char
 248 save(ImlibImage * im, ImlibProgressFunction progress, char progress_granularity)
 249 {
 250    FILE               *f;
 251    DATA32             *ptr;
 252    int                 y, pl = 0, alpha = 0;
 253    char                pper = 0;
 254
 255
 256    /* no image data? abort */
 257    if (!im->data)
 258       return 0;
 259    f = fopen(im->real_file, "wb");
 260    if (!f)
 261       return 0;
 262    if (im->flags & F_HAS_ALPHA)
 263       alpha = 1;
 264    fprintf(f, "ARGB %i %i %i\n", im->w, im->h, alpha);
 265    ptr = im->data;
 266    for (y = 0; y < im->h; y++)
 267      {
 268         fwrite(ptr, im->w, 4, f);
 269         ptr += im->w;
 270         if (progress)
 271           {
 272              char                per;
 273              int                 l;
 274
 275              per = (char)((100 * y) / im->h);
 276              if (((per - pper) >= progress_granularity) || (y == (im->h - 1)))
 277                {
 278                   l = y - pl;
 279                   if (!progress(im, per, 0, (y - l), im->w, l))
 280                     {
 281                        fclose(f);
 282                        return 2;
 283                     }
 284                   pper = per;
 285                   pl = y;
 286                }
 287           }
 288      }
 289    /* finish off */
 290    fclose(f);
 291    return 1;
 292 }
 293
 294 void
 295 formats(ImlibLoader * l)
 296 {
 297    char               *list_formats[] = { "fits", "fts", "fit" };
 298
 299    {
 300       int                 i;
 301
 302       l->num_formats = (sizeof(list_formats) / sizeof(char *));
 303       l->formats = malloc(sizeof(char *) * l->num_formats);
 304       for (i = 0; i < l->num_formats; i++)
 305          l->formats[i] = strdup(list_formats[i]);
 306    }
 307 }