git.oblomov.eu Git - zanmar/loader_fits/blob - loader_fits.c

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