IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 14, 2011, 2:58:01 PM (15 years ago)
Author:
watersc1
Message:

merge from czw_branch of logflux code

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psLib/src/fits/psFitsScale.c

    r27417 r30636  
    100100}
    101101
     102   
    102103// Determine appropriate BSCALE and BZERO for an image, mapping the standard deviation to the nominated number
    103104// of bits
    … …  
    212213    switch (options->scaling) {
    213214      case PS_FITS_SCALE_STDEV_POSITIVE:
     215      case PS_FITS_SCALE_LOG_STDEV_POSITIVE:
    214216        // Put (mean - N sigma) at the lowest possible value: predominantly positive images
    215217        imageVal = mean - options->stdevNum * stdev;
    … …  
    217219        break;
    218220      case PS_FITS_SCALE_STDEV_NEGATIVE:
     221      case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
    219222        // Put (mean + N sigma) at the highest possible value: predominantly negative images
    220223        imageVal = mean + options->stdevNum * stdev;
    … …  
    222225        break;
    223226      case PS_FITS_SCALE_STDEV_BOTH:
     227      case PS_FITS_SCALE_LOG_STDEV_BOTH:
    224228        // Put mean right in the middle: images with an equal abundance of positive and negative values
    225229        imageVal = mean;
    … …  
    237241}
    238242
     243
     244static bool logscaleStdev(double *bscale, // Scaling, to return
     245                          double *bzero, // Zero point, to return
     246                          double *boffset, // Log offset, to return
     247                          const psImage *image, // Image to scale
     248                          const psImage *mask, // Mask image
     249                          psImageMaskType maskVal, // Value to mask
     250                          const psFitsOptions *options // FITS options
     251    )
     252{
     253    psAssert(bscale, "impossible");
     254    psAssert(bzero, "impossible");
     255    psAssert(boffset, "impossible");
     256    psAssert(image, "impossible");
     257    psAssert(options, "impossible");
     258
     259    psTrace("psLib.fits", 3, "Scaling image by logarithm statistics");
     260    int numCols = image->numCols, numRows = image->numRows; // Size of image
     261
     262    psImage *copy;
     263   
     264    *boffset = 99e99;
     265
     266    // Make a copy of the image to pass to get the scaling parameters
     267
     268    switch (image->type.type) {
     269    case PS_TYPE_F32:
     270      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     271      break;
     272    case PS_TYPE_F64:
     273      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     274      break;
     275    default:
     276      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     277      return NULL;
     278      break;
     279    }
     280   
     281    // Determine the minimum value on this image.
     282    switch (image->type.type) {
     283    case PS_TYPE_F32:
     284      for (int y = 0; y < numRows; y++) {
     285        for (int x = 0; x < numCols; x++) {
     286          psF32 value = image->data.F32[y][x];
     287          if (isfinite(value)) {
     288            if (value < *boffset) {
     289              *boffset = value;
     290            }
     291          }
     292        }
     293      }
     294      break;
     295    case PS_TYPE_F64:
     296      for (int y = 0; y < numRows; y++) {
     297        for (int x = 0; x < numCols; x++) {
     298          psF64 value = image->data.F64[y][x];
     299          if (isfinite(value)) {
     300            if (value < *boffset) {
     301              *boffset = value;
     302            }
     303          }
     304        }
     305      }
     306      break;
     307    default:
     308      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     309      return NULL;
     310      break;
     311    }
     312    // We only need to offset images that go negative.
     313    if (*boffset > 0.0) {
     314      *boffset = 0.0;
     315    }
     316    // Write offset to header
     317    // How?
     318    //    psMetadataAddF32(header,PS_LIST_TAIL,"LOGZERO",0,"Flux offset subtracted before taking logarithm.",offset);
     319    // Take the logarithm of the image, applying the offset
     320    switch (image->type.type) {
     321    case PS_TYPE_F32:
     322      for (int y = 0; y < numRows; y++) {
     323        for (int x = 0; x < numCols; x++) {
     324/*        if (x == 2331 && y == 2843) { */
     325/*          fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */
     326/*        } */
     327          copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset));
     328        }
     329      }
     330      break;
     331    case PS_TYPE_F64:
     332        for (int y = 0; y < numRows; y++) {
     333          for (int x = 0; x < numCols; x++) {
     334            //      fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset));
     335            copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset));
     336          }
     337        }
     338      break;
     339    default:
     340      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     341      return NULL;
     342      break;
     343    }
     344     
     345    // Do regular scaling on the logarithm image
     346    if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) {
     347      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     348      return false;
     349    }
     350    psFree(copy);
     351    return true;
     352}
     353
     354static bool logscaleRange(double *bscale, // Scaling, to return
     355                          double *bzero, // Zero point, to return
     356                          double *boffset, // Log offset, to return
     357                          const psImage *image, // Image to scale
     358                          const psFitsOptions *options // FITS options
     359    )
     360{
     361    psAssert(bscale, "impossible");
     362    psAssert(bzero, "impossible");
     363    psAssert(boffset, "impossible");
     364    psAssert(image, "impossible");
     365    psAssert(options, "impossible");
     366
     367    psTrace("psLib.fits", 3, "Scaling image by logarithm statistics");
     368    int numCols = image->numCols, numRows = image->numRows; // Size of image
     369
     370    psImage *copy;
     371   
     372    *boffset = 99e99;
     373
     374    // Make a copy of the image to pass to get the scaling parameters
     375
     376    switch (image->type.type) {
     377    case PS_TYPE_F32:
     378      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     379      break;
     380    case PS_TYPE_F64:
     381      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     382      break;
     383    default:
     384      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     385      return NULL;
     386      break;
     387    }
     388   
     389    // Determine the minimum value on this image.
     390    switch (image->type.type) {
     391    case PS_TYPE_F32:
     392      for (int y = 0; y < numRows; y++) {
     393        for (int x = 0; x < numCols; x++) {
     394          psF32 value = image->data.F32[y][x];
     395          if (!isfinite(value)) {
     396            if (value < *boffset) {
     397              *boffset = value;
     398            }
     399          }
     400        }
     401      }
     402      break;
     403    case PS_TYPE_F64:
     404      for (int y = 0; y < numRows; y++) {
     405        for (int x = 0; x < numCols; x++) {
     406          psF64 value = image->data.F64[y][x];
     407          if (!isfinite(value)) {
     408            if (value < *boffset) {
     409              *boffset = value;
     410            }
     411          }
     412        }
     413      }
     414      break;
     415    default:
     416      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     417      return NULL;
     418      break;
     419    }
     420    // We only need to offset images that go negative.
     421    if (*boffset > 0.0) {
     422      *boffset = 0.0;
     423    }
     424    // Write offset to header
     425    // How?
     426    //    psMetadataAddF32(header,PS_LIST_TAIL,"LOGZERO",0,"Flux offset subtracted before taking logarithm.",offset);
     427    // Take the logarithm of the image, applying the offset
     428    switch (image->type.type) {
     429    case PS_TYPE_F32:
     430      for (int y = 0; y < numRows; y++) {
     431        for (int x = 0; x < numCols; x++) {
     432          copy->data.F32[y][x] = (log10( image->data.F32[y][x] - *boffset));
     433        }
     434      }
     435      break;
     436    case PS_TYPE_F64:
     437        for (int y = 0; y < numRows; y++) {
     438          for (int x = 0; x < numCols; x++) {
     439            copy->data.F64[y][x] = (log10( image->data.F64[y][x] - *boffset));
     440          }
     441        }
     442      break;
     443    default:
     444      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     445      return NULL;
     446      break;
     447    }
     448     
     449    // Do regular scaling on the logarithm image
     450    if (!scaleRange(bscale, bzero, copy, options)) {
     451      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     452      return false;
     453    }
     454    psFree(copy);
     455    return true;
     456}
     457
     458
    239459//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    240460// Public functions
    241461//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    242462
    243 bool psFitsScaleDetermine(double *bscale, double *bzero, long *blank, const psImage *image,
     463bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image,
    244464                          const psImage *mask, psImageMaskType maskVal, const psFits *fits)
    245465{
    246466    PS_ASSERT_PTR_NON_NULL(bscale, false);
    247467    PS_ASSERT_PTR_NON_NULL(bzero, false);
     468    PS_ASSERT_PTR_NON_NULL(boffset, false);
    248469    PS_ASSERT_PTR_NON_NULL(blank, false);
    249470    PS_ASSERT_IMAGE_NON_NULL(image, false);
    … …  
    256477    *bscale = NAN;
    257478    *bzero = NAN;
     479    *boffset = 0;
    258480    *blank = 0;
    259481
    … …  
    299521        }
    300522        break;
     523    case PS_FITS_SCALE_LOG_RANGE:
     524      if (!logscaleRange(bscale,bzero,boffset,image,options)) {
     525        psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range");
     526        return false;
     527      }
     528      break;
     529    case PS_FITS_SCALE_LOG_STDEV_POSITIVE:
     530    case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
     531    case PS_FITS_SCALE_LOG_STDEV_BOTH:
     532      if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) {
     533            psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     534            return false;
     535        }
     536        break;
    301537      case PS_FITS_SCALE_MANUAL:
    302538        *bscale = options->bscale;
    303539        *bzero = options->bzero;
     540        break;
     541      case PS_FITS_SCALE_LOG_MANUAL:
     542        *bscale = options->bscale;
     543        *bzero = options->bzero;
     544        *boffset = options->boffset;
    304545        break;
    305546      default:
    … …  
    316557
    317558
    318     psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BLANK = %ld\n", *bscale, *bzero, *blank);
     559    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n",
     560            *bscale, *bzero, *boffset, *blank);
    319561    return true;
    320562}
    321563
    322564
    323 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero,
     565psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset,
    324566                            psRandom *rng)
    325567{
    … …  
    377619        for (int y = 0; y < numRows; y++) { \
    378620            for (int x = 0; x < numCols; x++) { \
    379                 ps##INTYPE value = (IN)->data.INTYPE[y][x]; \
     621              ps##INTYPE value; \
     622              if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||       \
     623                (options->scaling == PS_FITS_SCALE_LOG_MANUAL)||        \
     624                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
     625                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
     626                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
     627                value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
     628              }                                                         \
     629                else { \
     630                  value = (IN)->data.INTYPE[y][x];      \
     631                }                                           \
    380632                if (!isfinite(value)) { \
    381633                    /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
    … …  
    423675
    424676
    425 #if 0
     677
    426678// This function to apply BSCALE and BZERO to an image read immediately from disk should not be necessary at
    427679// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
    428680// later desire it (e.g., if we ever make our own FITS implementation).
    429 psImage *psFitsScaleFromDisk(psFits *fits, psImage *image)
     681psImage *psFitsScaleFromDisk(const psImage *image, double boffset)
    430682{
    431683    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    432 
    433     if (bscale == 0.0) {
    434         // BSCALE = 0 means don't apply anything
    435         return psMemIncrRefCounter(image);
    436     }
    437684
    438685    psElemType inType = image->type.type; // Type for input image
    … …  
    471718      for (int y = 0; y < numRows; y++) { \
    472719          for (int x = 0; x < numCols; x++) { \
    473               out->data.OUTTYPE[y][x] = image->data.INTYPE[y][x] * bscale + bzero; \
     720            out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \
    474721          } \
    475722      } \
    … …  
    505752    return out;
    506753}
    507 #endif
    508 
    509 
    510754
    511755psFitsScaling psFitsScalingFromString(const char *string)
    … …  
    517761    if (strcasecmp(string, "STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_STDEV_NEGATIVE;
    518762    if (strcasecmp(string, "STDEV_BOTH") == 0)     return PS_FITS_SCALE_STDEV_BOTH;
     763    if (strcasecmp(string, "LOG_RANGE") == 0)      return PS_FITS_SCALE_LOG_RANGE;
     764    if (strcasecmp(string, "LOG_MANUAL") == 0)      return PS_FITS_SCALE_LOG_MANUAL;
     765    if (strcasecmp(string, "LOG_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_POSITIVE;
     766    if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE;
     767    if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH;
    519768    if (strcasecmp(string, "MANUAL") == 0)         return PS_FITS_SCALE_MANUAL;
    520769
Note: See TracChangeset for help on using the changeset viewer.