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