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 void getRGB( double, double, double, short, short, DATA8 *, DATA8 *, DATA8 * );
  17
  18 char
  19 load(ImlibImage * im, ImlibProgressFunction progress,
  20      char progress_granularity, char immediate_load)
  21 {
  22
  23    int                 w = 0, h = 0, alpha = 0;
  24    //FILE               *f;
  25    fitsfile            *f;
  26    int                 status = 0;
  27
  28    if (im->data)
  29       return 0;
  30    //f = fopen(im->real_file, "rb");
  31    if ( fits_open_file( &f, im->real_file, READONLY, &status) ) {
  32          return 0;
  33    }
  34
  35    //if (!f)
  36    //   return 0;
  37
  38    /* header */
  39    {
  40       long                naxes[2];
  41       int                 nfound;
  42       //char                buf[256], buf2[256];
  43
  44       /* read the NAXIS1 and NAXIS2 keyword to get image size */
  45       if ( fits_read_keys_lng( f, "NAXIS", 1, 2, naxes, &nfound, &status ) ){
  46            fits_close_file(f, &status);
  47            return 0;
  48       }
  49       w = naxes[ 0 ];
  50       h = naxes[ 1 ];
  51
  52       /*
  53       buf[0] = '\0';
  54       if (!fgets(buf, 255, f))
  55         {
  56            fclose(f);
  57            return 0;
  58         }
  59       buf2[0] = '\0';
  60       sscanf(buf, "%s %i %i %i", buf2, &w, &h, &alpha);
  61       if (strcmp(buf2, "ARGB"))
  62         {
  63            fclose(f);
  64            return 0;
  65         }
  66         */
  67       if (!IMAGE_DIMENSIONS_OK(w, h))
  68         {
  69            fits_close_file(f, &status);
  70            //fclose(f);
  71            return 0;
  72         }
  73       im->w = w;
  74       im->h = h;
  75       if (!im->format)
  76         {
  77            if (alpha)
  78               SET_FLAG(im->flags, F_HAS_ALPHA);
  79            else
  80               UNSET_FLAG(im->flags, F_HAS_ALPHA);
  81            im->format = strdup("fits");
  82         }
  83    }
  84    if (((!im->data) && (im->loader)) || (immediate_load) || (progress))
  85      {
  86         DATA32             *ptr;
  87         DATA8              r, g, b;
  88         int                y, x, pl = 0;
  89         char               pper = 0;
  90         double             *buffer, datamin, datamax, nullval = 0;
  91         int                anynull, nbuffer;
  92         long               fpixel = 1;
  93         long               *bin, nmax = 0, ndataok, sumpix = 0;
  94         int                nbin, k, kmin, mostfreq = 0;
  95         double             binw;
  96         char               *envdmin, *envdmax;
  97         short              color, inv;
  98
  99     fprintf(stderr, "%s\n", im->real_file );
 100         /* must set the im->data member before callign progress function */
 101         ptr = im->data = malloc(w * h * sizeof(DATA32));
 102         nbuffer = w;
 103         buffer = malloc( nbuffer * sizeof( double ) );
 104
 105         if (!im->data)
 106           {
 107              fits_close_file(f, &status);
 108              //fclose(f);
 109              return 0;
 110           }
 111
 112         envdmin = getenv( "DATAMIN" );
 113         envdmax = getenv( "DATAMAX" );
 114         if( envdmin != NULL && envdmax != NULL ){
 115             datamin = atof( envdmin );
 116             datamax = atof( envdmax );
 117             fprintf(stderr, "ENV:data min %g, max %g\n", datamin, datamax);}
 118         else{
 119
 120             datamin = HUGE_VALF;
 121             datamax = -HUGE_VALF;
 122
 123             //Find min max
 124             for (y = 0; y < h; y++)
 125               {
 126                  fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 127                               buffer, &anynull, &status);
 128                  for( x = 0; x < w; x++ ){
 129                      if( buffer[ x ] > datamax )
 130                          datamax = buffer[ x ];
 131                      if( buffer[ x ] < datamin )
 132                          datamin = buffer[ x ];
 133                  }
 134                  fpixel += nbuffer;
 135               }
 136             fprintf(stderr, "data min %g, max %g\n", datamin, datamax);
 137
 138             nbin = 4 * sqrt( datamax - datamin );
 139             if( nbin > 5000 ) nbin = 1000;
 140             binw = ( datamax - datamin ) / nbin;
 141             fprintf(stderr, "nbin and bin width:%d %g\n", nbin, binw);
 142             bin = calloc( nbin, sizeof( long ) );
 143
 144             //Find histogram
 145             fpixel = 1;
 146             for (y = 0; y < h; y++)
 147               {
 148                  fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 149                               buffer, &anynull, &status);
 150                  for( x = 0; x < w; x++ ){
 151                      for( k = 0; k < nbin; k++ ){
 152                          if( buffer[ x ] >= k * binw &&
 153                              buffer[ x ] < ( k + 1 ) * binw ){
 154                              bin[ k ] += 1;
 155                              break;
 156                          }
 157                      }
 158                  }
 159                  fpixel += nbuffer;
 160               }
 161             //Find most frequent pixel value
 162             for( k = 0; k < nbin; k++ ){
 163                 if( bin[ k ] > nmax ){
 164                     nmax = bin[ k ];
 165                     mostfreq = k;
 166                 }
 167                 //fprintf(stderr, "%ld ", bin[ k ] );
 168             }
 169             fprintf(stderr, "\nmostfreq and nmax = %d %ld ", mostfreq, nmax );
 170
 171             //Set new DATAMIN. look for FWHM on a log scale
 172             nmax = log( nmax ) / 2;
 173             datamin = ( mostfreq + 0.5 ) * binw;
 174             for( k = mostfreq - 1; k >= 0; k-- ){
 175                 if( bin[ k ] > 0 )
 176                     if( log( bin[ k ] ) < nmax ){
 177                         kmin = k + 1;
 178                         datamin = ( (kmin) + 0.5 ) * binw;
 179                         break;
 180                     }
 181             }
 182             if( k == 0 ){   //case were we have one single peak surrounded by 0
 183                kmin = mostfreq;
 184                datamin = ( (kmin) + 0.5 ) * binw;
 185             }
 186
 187             //Set new DATAMAX when we have 99.5% of the pixels accounted for
 188             ndataok = h * w;
 189             for( k = 0; k < kmin; k++)
 190                 ndataok -= bin[ k ];
 191             //ndataok = h * w; //To get values closer to ds9's uncomment this
 192             ndataok = 0.997 * ndataok;
 193             for( k = kmin; k < nbin; k++){
 194                 sumpix += bin[ k ];
 195                 if( sumpix > ndataok )
 196                     break;
 197             }
 198             datamax = ( ( k + 1 ) + 0.5 ) * binw;
 199         }
 200
 201         free( bin );
 202
 203         // use the new datamin datamax values to rescale and create imlib im
 204         ptr += w * h - 1;
 205         fpixel = 1;
 206         color = 1;
 207         inv = 0;
 208         for (y = 0; y < h; y++)
 209           {
 210              fits_read_img( f, TDOUBLE, fpixel, nbuffer, &nullval,
 211                           buffer, &anynull, &status);
 212
 213              for( x = w - 1; x >= 0; x-- ){
 214                  getRGB( buffer[ x ], datamin, datamax, color, inv, &r, &g, &b);
 215                  if( x == w / 100 && y == h / 100 ) fprintf( stderr, "\n%d %d %d %g\n",r,g,b, buffer[ x ]);
 216                  *ptr = (((( (0xff << 8) | r) << 8) | g) << 8) | b ;
 217                  ptr--;
 218              }
 219              fpixel += nbuffer;
 220
 221              if (progress)
 222                {
 223                   char                per;
 224                   int                 l;
 225
 226                   per = (char)((100 * y) / im->h);
 227                   if (((per - pper) >= progress_granularity) ||
 228                       (y == (im->h - 1)))
 229                     {
 230                        l = y - pl;
 231                        if (!progress(im, per, 0, (y - l), im->w, l))
 232                          {
 233                             fits_close_file(f, &status);
 234                             //fclose(f);
 235                             return 2;
 236                          }
 237                        pper = per;
 238                        pl = y;
 239                     }
 240                }
 241           }
 242           fprintf(stderr, "\ndata min %g, max %g\n", datamin, datamax);
 243           free( buffer );
 244           //free( ptr );
 245      }
 246    //fclose(f);
 247    fits_close_file(f, &status);
 248    return 1;
 249 }
 250
 251 char
 252 save(ImlibImage * im, ImlibProgressFunction progress, char progress_granularity)
 253 {
 254    FILE               *f;
 255    DATA32             *ptr;
 256    int                 y, pl = 0, alpha = 0;
 257    char                pper = 0;
 258
 259
 260    /* no image data? abort */
 261    if (!im->data)
 262       return 0;
 263    f = fopen(im->real_file, "wb");
 264    if (!f)
 265       return 0;
 266    if (im->flags & F_HAS_ALPHA)
 267       alpha = 1;
 268    fprintf(f, "ARGB %i %i %i\n", im->w, im->h, alpha);
 269    ptr = im->data;
 270    for (y = 0; y < im->h; y++)
 271      {
 272         fwrite(ptr, im->w, 4, f);
 273         ptr += im->w;
 274         if (progress)
 275           {
 276              char                per;
 277              int                 l;
 278
 279              per = (char)((100 * y) / im->h);
 280              if (((per - pper) >= progress_granularity) || (y == (im->h - 1)))
 281                {
 282                   l = y - pl;
 283                   if (!progress(im, per, 0, (y - l), im->w, l))
 284                     {
 285                        fclose(f);
 286                        return 2;
 287                     }
 288                   pper = per;
 289                   pl = y;
 290                }
 291           }
 292      }
 293    /* finish off */
 294    fclose(f);
 295    return 1;
 296 }
 297
 298 void
 299 formats(ImlibLoader * l)
 300 {
 301    char               *list_formats[] = { "fits", "fts", "fit" };
 302
 303    {
 304       int                 i;
 305
 306       l->num_formats = (sizeof(list_formats) / sizeof(char *));
 307       l->formats = malloc(sizeof(char *) * l->num_formats);
 308       for (i = 0; i < l->num_formats; i++)
 309          l->formats[i] = strdup(list_formats[i]);
 310    }
 311 }
 312
 313 void getRGB( double val, double min, double max, short color, short inv,
 314         DATA8 *r, DATA8 *g, DATA8 *b ){
 315     double fac;
 316     short ind = 0;
 317     const float slsR[200]={0, 0.043442, 0.086883, 0.130325, 0.173767, 0.217208,
 318         0.260650, 0.304092, 0.347533, 0.390975, 0.434417, 0.477858, 0.521300,
 319         0.506742, 0.492183, 0.477625, 0.463067, 0.448508, 0.433950, 0.419392,
 320         0.404833, 0.390275, 0.375717, 0.361158, 0.346600, 0.317717, 0.288833,
 321         0.259950, 0.231067, 0.202183, 0.173300, 0.144417, 0.115533, 0.086650,
 322         0.057767, 0.028883, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 323         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 324         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.033167,
 325         0.066333, 0.099500, 0.132667, 0.165833, 0.199000, 0.232167, 0.265333,
 326         0.298500, 0.331667, 0.364833, 0.398000, 0.430950, 0.463900, 0.496850,
 327         0.529800, 0.562750, 0.595700, 0.628650, 0.661600, 0.694550, 0.727500,
 328         0.760450, 0.793400, 0.810617, 0.827833, 0.845050, 0.862267, 0.879483,
 329         0.896700, 0.913917, 0.931133, 0.948350, 0.965567, 0.982783, 1, 0.995725,
 330         0.991450, 0.987175, 0.982900, 0.978625, 0.974350, 0.970075, 0.965800,
 331         0.961525, 0.957250, 0.952975, 0.948700, 0.952975, 0.957250, 0.961525,
 332         0.965800, 0.970075, 0.974350, 0.978625, 0.982900, 0.987175, 0.991450,
 333         0.995725, 1, 0.998342, 0.996683, 0.995025, 0.993367, 0.991708, 0.990050,
 334         0.988392, 0.986733, 0.985075, 0.983417, 0.981758, 0.980100, 0.955925,
 335         0.931750, 0.907575, 0.883400, 0.859225, 0.835050, 0.810875, 0.786700,
 336         0.762525, 0.738350, 0.714175, 0.690000, 0.715833, 0.741667, 0.767500,
 337         0.793333, 0.819167, 0.845000, 0.870833, 0.896667, 0.922500, 0.948333,
 338         0.974167, 1, 1, 1, 1, 1, 1, 1, 1};
 339
 340     const float slsG[200]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 341         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.019817,
 342         0.039633, 0.059450, 0.079267, 0.099083, 0.118900, 0.138717, 0.158533,
 343         0.178350, 0.198167, 0.217983, 0.237800, 0.268533, 0.299267, 0.330000,
 344         0.360733, 0.391467, 0.422200, 0.452933, 0.483667, 0.514400, 0.545133,
 345         0.575867, 0.606600, 0.631733, 0.656867, 0.682000, 0.707133, 0.732267,
 346         0.757400, 0.782533, 0.807667, 0.832800, 0.857933, 0.883067, 0.908200,
 347         0.901908, 0.895617, 0.889325, 0.883033, 0.876742, 0.870450, 0.864158,
 348         0.857867, 0.851575, 0.845283, 0.838992, 0.832700, 0.832308, 0.831917,
 349         0.831525, 0.831133, 0.830742, 0.830350, 0.829958, 0.829567, 0.829175,
 350         0.828783, 0.828392, 0.828000, 0.834167, 0.840333, 0.846500, 0.852667,
 351         0.858833, 0.865000, 0.871167, 0.877333, 0.883500, 0.889667, 0.895833,
 352         0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000,
 353         0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.902000, 0.897133,
 354         0.892267, 0.887400, 0.882533, 0.877667, 0.872800, 0.867933, 0.863067,
 355         0.858200, 0.853333, 0.848467, 0.843600, 0.824892, 0.806183, 0.787475,
 356         0.768767, 0.750058, 0.731350, 0.712642, 0.693933, 0.675225, 0.656517,
 357         0.637808, 0.619100, 0.600408, 0.581717, 0.563025, 0.544333, 0.525642,
 358         0.506950, 0.488258, 0.469567, 0.450875, 0.432183, 0.413492, 0.394800,
 359         0.361900, 0.329000, 0.296100, 0.263200, 0.230300, 0.197400, 0.164500,
 360         0.131600, 0.098700, 0.065800, 0.032900, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 361         0, 0, 0, 0.083333, 0.166667, 0.250000, 0.333333, 0.416667, 0.500000,
 362         0.583333, 0.666667, 0.750000, 0.833333, 0.916667, 1, 1, 1, 1, 1, 1, 1,
 363         1};
 364
 365     const float slsB[200]={0, 0.052883, 0.105767, 0.158650, 0.211533, 0.264417,
 366         0.317300, 0.370183, 0.423067, 0.475950, 0.528833, 0.581717, 0.634600,
 367         0.640217, 0.645833, 0.651450, 0.657067, 0.662683, 0.668300, 0.673917,
 368         0.679533, 0.685150, 0.690767, 0.696383, 0.702000, 0.712192, 0.722383,
 369         0.732575, 0.742767, 0.752958, 0.763150, 0.773342, 0.783533, 0.793725,
 370         0.803917, 0.814108, 0.824300, 0.838942, 0.853583, 0.868225, 0.882867,
 371         0.897508, 0.912150, 0.926792, 0.941433, 0.956075, 0.970717, 0.985358,
 372         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.975300, 0.950600, 0.925900,
 373         0.901200, 0.876500, 0.851800, 0.827100, 0.802400, 0.777700, 0.753000,
 374         0.728300, 0.703600, 0.676675, 0.649750, 0.622825, 0.595900, 0.568975,
 375         0.542050, 0.515125, 0.488200, 0.461275, 0.434350, 0.407425, 0.380500,
 376         0.354858, 0.329217, 0.303575, 0.277933, 0.252292, 0.226650, 0.201008,
 377         0.175367, 0.149725, 0.124083, 0.098442, 0.072800, 0.066733, 0.060667,
 378         0.054600, 0.048533, 0.042467, 0.036400, 0.030333, 0.024267, 0.018200,
 379         0.012133, 0.006067, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.003983,
 380         0.007967, 0.011950, 0.015933, 0.019917, 0.023900, 0.027883, 0.031867,
 381         0.035850, 0.039833, 0.043817, 0.047800, 0.051600, 0.055400, 0.059200,
 382         0.063000, 0.066800, 0.070600, 0.074400, 0.078200, 0.082000, 0.085800,
 383         0.089600, 0.093400, 0.085617, 0.077833, 0.070050, 0.062267, 0.054483,
 384         0.046700, 0.038917, 0.031133, 0.023350, 0.015567, 0.007783, 0, 0, 0,
 385         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 386         0.083333, 0.166667, 0.250000, 0.333333, 0.416667, 0.500000, 0.583333,
 387         0.666667, 0.750000, 0.833333, 0.916667, 1, 1, 1, 1, 1, 1, 1, 1};
 388
 389     if( val <= min )
 390         val = min;
 391     else
 392         if( val >= max )
 393             val = max;
 394
 395     fac = ( val - min ) / ( max - min );
 396
 397     if( inv ) fac = 1 - fac;
 398
 399     if( color == 0 ){ // gray
 400         *r = ceil ( fac * 255 );
 401         *g = *r;
 402         *b = *r;
 403         return;
 404     }
 405     if( color == 1 ){ // color SLS
 406         ind = ceil( fac * 199 );    // color table has 200 elements
 407         *r = ceil( slsR[ ind ] * 255 );
 408         *g = ceil( slsG[ ind ] * 255 );
 409         *b = ceil( slsB[ ind ] * 255 );
 410         return;
 411     }
 412
 413
 414         return;
 415 }
 416
 417