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