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( log( bin[ k ] ) < nmax ){
 173                     kmin = k + 1;
 174                     datamin = ( (kmin) + 0.5 ) * binw;
 175                     break;
 176                 }
 177             }
 178             //Set new DATAMAX when we have 99.5% of the pixels accounted for
 179             ndataok = h * w;
 180             for( k = 0; k < kmin; k++)
 181                 ndataok -= bin[ k ];
 182             //ndataok = h * w; //To get values closer to ds9's uncomment this
 183             for( k = kmin; k < nbin; k++){
 184                 sumpix += bin[ k ];
 185                 if( 1.0 * sumpix / ndataok > 0.997 )
 186                     break;
 187             }
 188             datamax = ( ( k + 1 ) + 0.5 ) * binw;
 189         }
 190
 191         // use the new datamin datamax values to rescale and create imlib im
 192         ptr += w * h - 1;
 193         fpixel = 1;
 194         for (y = 0; y < h; y++)
 195           {
 196              fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 197                           buffer, &anynull, &status);
 198
 199              for( x = w - 1; x >= 0; x-- ){
 200                  if( buffer[ x ] <= datamin )
 201                      gray = 0;
 202                  else
 203                      if( buffer[ x ] >= datamax )
 204                          gray = 255;
 205                      else
 206                          gray = ceil(( buffer[ x ] - datamin ) / ( datamax - datamin ) * 255 );
 207                  *ptr = (((( (0xff << 8) | gray) << 8) | gray) << 8) | gray ;
 208                  ptr--;
 209              }
 210              fpixel += nbuffer;
 211
 212              if (progress)
 213                {
 214                   char                per;
 215                   int                 l;
 216
 217                   per = (char)((100 * y) / im->h);
 218                   if (((per - pper) >= progress_granularity) ||
 219                       (y == (im->h - 1)))
 220                     {
 221                        l = y - pl;
 222                        if (!progress(im, per, 0, (y - l), im->w, l))
 223                          {
 224                             fits_close_file(f, &status);
 225                             //fclose(f);
 226                             return 2;
 227                          }
 228                        pper = per;
 229                        pl = y;
 230                     }
 231                }
 232           }
 233             fprintf(stderr, "data min %g, max %g\n", datamin, datamax);
 234      }
 235    //fclose(f);
 236    fits_close_file(f, &status);
 237    return 1;
 238 }
 239
 240 char
 241 save(ImlibImage * im, ImlibProgressFunction progress, char progress_granularity)
 242 {
 243    FILE               *f;
 244    DATA32             *ptr;
 245    int                 y, pl = 0, alpha = 0;
 246    char                pper = 0;
 247
 248
 249    /* no image data? abort */
 250    if (!im->data)
 251       return 0;
 252    f = fopen(im->real_file, "wb");
 253    if (!f)
 254       return 0;
 255    if (im->flags & F_HAS_ALPHA)
 256       alpha = 1;
 257    fprintf(f, "ARGB %i %i %i\n", im->w, im->h, alpha);
 258    ptr = im->data;
 259    for (y = 0; y < im->h; y++)
 260      {
 261         fwrite(ptr, im->w, 4, f);
 262         ptr += im->w;
 263         if (progress)
 264           {
 265              char                per;
 266              int                 l;
 267
 268              per = (char)((100 * y) / im->h);
 269              if (((per - pper) >= progress_granularity) || (y == (im->h - 1)))
 270                {
 271                   l = y - pl;
 272                   if (!progress(im, per, 0, (y - l), im->w, l))
 273                     {
 274                        fclose(f);
 275                        return 2;
 276                     }
 277                   pper = per;
 278                   pl = y;
 279                }
 280           }
 281      }
 282    /* finish off */
 283    fclose(f);
 284    return 1;
 285 }
 286
 287 void
 288 formats(ImlibLoader * l)
 289 {
 290    char               *list_formats[] = { "fits", "fts", "fit" };
 291
 292    {
 293       int                 i;
 294
 295       l->num_formats = (sizeof(list_formats) / sizeof(char *));
 296       l->formats = malloc(sizeof(char *) * l->num_formats);
 297       for (i = 0; i < l->num_formats; i++)
 298          l->formats[i] = strdup(list_formats[i]);
 299    }
 300 }