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