IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31618


Ignore:
Timestamp:
Jun 8, 2011, 2:41:13 PM (15 years ago)
Author:
watersc1
Message:

compilable version of updated stack code (doesn't correctly reject) and asinh scaling (completely untested).

Location:
branches/czw_branch/20110406
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20110406/ippconfig/recipes/fitstypes.mdc

    r30636 r31618  
    101101COMP_STACK     METADATA
    102102        BITPIX          S32     16
    103         SCALING         STR     LOG_STDEV_POSITIVE
     103        SCALING         STR     ASINH_STDEV_POSITIVE
    104104        STDEV.BITS      S32     4
    105105        STDEV.NUM       F32     10

  • branches/czw_branch/20110406/ippconfig/recipes/reductionClasses.mdc

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/czw_branch/20110406/psLib/src/fits/psFits.h

    r30636 r31618  
    5555    PS_FITS_SCALE_LOG_STDEV_NEGATIVE,   ///< Take logarithm, Auto-scale to sample stdev, place mean at upper limit
    5656    PS_FITS_SCALE_LOG_STDEV_BOTH,       ///< Take logarithm, Auto-scale to sample stdev, place mean at middle
    57     PS_FITS_SCALE_LOG_MANUAL            ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET)
     57    PS_FITS_SCALE_LOG_MANUAL,           ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET)
     58    PS_FITS_SCALE_ASINH_RANGE,          ///< Do asinh scaling, auto-scale to preserve dynamic range
     59    PS_FITS_SCALE_ASINH_STDEV_POSITIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at lower limit
     60    PS_FITS_SCALE_ASINH_STDEV_NEGATIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at upper limit
     61    PS_FITS_SCALE_ASINH_STDEV_BOTH,     ///< Do asinh scaling, auto-scale to sample stdev, place mean at middle
     62    PS_FITS_SCALE_ASINH_MANUAL          ///< Manual scaling after doing asinh scaling.
    5863} psFitsScaling;
    5964
    … …  
    7277    double bscale, bzero;               ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL)
    7378    double boffset;                     ///< Manually specified BOFFSET (for SCALE_MANUAL)
     79    double bsoften;                     ///< Manuallly specified BSOFTEN (for SCALE_MANUAL)
    7480    double mean, stdev;                 ///< Mean and standard deviation of image
    7581    int stdevBits;                      ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*)

  • branches/czw_branch/20110406/psLib/src/fits/psFitsImage.c

    r31606 r31618  
    5555    int psDatatype;                     // psLib data type
    5656    bool is_logscaled;                  // is this image log scaled using BOFFSET?
     57    bool is_asinh;                      // is this image asinh scaled using BSOFTEN?
    5758} p_psFitsReadInfo;
    5859
    … …  
    129130
    130131    // Check scale and zero
    131     double bscale = 0.0, bzero = 0.0, boffset = NAN;    // Scale and zero point
     132    double bscale = 0.0, bzero = 0.0, boffset = NAN, bsoften = NAN;    // Scale and zero point
    132133    if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) {
    133134        psFitsError(status, true, "Unable to read header.");
    … …  
    153154      info->is_logscaled = false;
    154155    }
     156    status = 0;
     157    if (fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status) && status != KEY_NO_EXIST) {
     158        psFitsError(status, true, "Unable to read header.");
     159        goto bad;
     160    }
     161    if (status == KEY_NO_EXIST) {
     162      info->is_asinh = false;
     163    }
     164    else if (isfinite(bsoften)) {
     165      info->is_asinh = true;
     166    }
     167    else {
     168      info->is_asinh = false;
     169    }
     170
    155171    status = 0;
    156172
    … …  
    262278                                          double *bzero, // Zero point applied
    263279                                          double *boffset, // Log offset applied
     280                                          double *bsoften, // asinh offset applied
    264281                                          long *blank, // Blank value (integer data)
    265282                                          psFitsFloat *floatType, // Type of custom floating-point
    … …  
    275292    psAssert(bzero, "impossible");
    276293    psAssert(boffset, "impossible");
     294    psAssert(bsoften, "impossible");
    277295    psAssert(floatType, "impossible");
    278296    psAssert(fits, "impossible");
    … …  
    322340        if (newScaleZero) {
    323341            // Choose an appropriate BSCALE and BZERO
    324           if (!psFitsScaleDetermine(bscale, bzero, boffset, blank, image, mask, maskVal, fits)) {
     342          if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) {
    325343                // We can't have the write dying for this reason --- try to save it somehow!
    326344                psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
    … …  
    348366        }
    349367
    350         return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng);
     368        return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng);
    351369    }
    352370
    … …  
    487505      status = 0;
    488506      fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status);
    489       psImage *newImage = psFitsScaleFromDisk(outImage,boffset);
     507      psImage *newImage = psFitsScaleFromDisk(outImage,boffset,0.0);
     508      psFree (outImage);
     509      outImage = newImage;
     510    }
     511    else if (info->is_asinh) {
     512      double bsoften;
     513      int status;
     514      status = 0;
     515      fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status);
     516      psImage *newImage = psFitsScaleFromDisk(outImage,0.0,bsoften);
    490517      psFree (outImage);
    491518      outImage = newImage;
    … …  
    602629    double bscale = NAN, bzero = NAN;   // Scale and zero point to put in header
    603630    double boffset = NAN;               // Log offset to put into header.
     631    double bsoften = NAN;               // Asinh softening parameter to put into header
    604632    long blank = 0;                     // Blank (undefined) value for image
    605633    psFitsFloat floatType;              // Custom floating-point convention type
    606     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, image,
     634    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image,
    607635                                                   mask, maskVal, NULL, true); // Image to write out
    608636    if (!diskImage) {
    … …  
    831859    double bscale = NAN, bzero = NAN;   // Scale and zero point to put in header
    832860    double boffset = NAN;               // Log offset to put in header
     861    double bsoften = NAN;               // Asinh softening parameter to put in header
    833862    long blank = 0;                     // Blank (undefined) value for image
    834863    psFitsFloat floatType;              // Custom floating-point convention type
    835     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, input,
     864    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input,
    836865                                                   mask, maskVal, NULL, false); // Image to write out
    837866    if (!diskImage) {

  • branches/czw_branch/20110406/psLib/src/fits/psFitsScale.c

    r30636 r31618  
    242242
    243243
     244   
     245
    244246static bool logscaleStdev(double *bscale, // Scaling, to return
    245247                          double *bzero, // Zero point, to return
    … …  
    456458}
    457459
     460static bool asinhStdev(double *bscale, // Scaling, to return
     461                       double *bzero,  // Zero point, to return
     462                       double *bsoften, // asinh softening parameter, to return
     463                       const psImage *image, // Image to scale
     464                       const psImage *mask,  // Mask image
     465                       psImageMaskType maskVal, // value to mask
     466                       const psFitsOptions *options // FITS options
     467                       )
     468{
     469  psAssert(bscale, "impossible");
     470  psAssert(bzero, "impossible");
     471  psAssert(bsoften, "impossible");
     472  psAssert(image, "impossible");
     473  psAssert(options, "impossible");
     474
     475  psTrace("psLib.fits", 3, "Scaling image by asinh method");
     476  int numCols = image->numCols, numRows = image->numRows; // Size of image
     477 
     478    // Measure the mean and stdev
     479    // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a
     480    // mask.
     481    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     482    psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object
     483    double mean, stdev;                                    // Mean and standard deviation
     484    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
     485        // It could be because the image is entirely masked, in which case we don't want to error
     486        bool good = false;              // Any good pixels?
     487
     488
     489// Find good pixels in an image, by image type
     490#define GOOD_PIXELS_CASE(TYPE) \
     491      case PS_TYPE_##TYPE: \
     492        for (int y = 0; y < image->numRows && !good; y++) { \
     493            for (int x = 0; x < image->numCols && !good; x++) { \
     494                if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \
     495                    continue; \
     496                } \
     497                if (!isfinite(image->data.TYPE[y][x])) { \
     498                    continue; \
     499                } \
     500                good = true; \
     501            } \
     502        } \
     503        break;
     504
     505        switch (image->type.type) {
     506            GOOD_PIXELS_CASE(F32);
     507            GOOD_PIXELS_CASE(F64);
     508          default:
     509            psAbort("Unsupported case: %x", image->type.type);
     510        }
     511
     512        if (!good) {
     513            psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0");
     514            psErrorClear();
     515            *bscale = 1.0;
     516            *bzero = 0.0;
     517            *bsoften = NAN;
     518            psFree(rng);
     519            psFree(stats);
     520            return true;
     521        }
     522
     523        // There are some good pixels in there somewhere; psImageBackground just didn't find them
     524        psLogMsg("psLib.fits", PS_LOG_DETAIL,
     525                 "Couldn't measure background statistics for image quantisation; retrying.");
     526        psErrorClear();
     527        // Retry using all the available pixels
     528        stats->nSubsample = image->numCols * image->numRows + 1;
     529        if (!psImageStats(stats, image, mask, maskVal)) {
     530            psLogMsg("psLib.fits", PS_LOG_DETAIL,
     531                     "Couldn't measure background statistics for image quantisation (attempt 2); retrying.");
     532            psErrorClear();
     533            // Retry with desperate statistic
     534            stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT;
     535            if (!psImageStats(stats, image, mask, maskVal)) {
     536                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image");
     537                psFree(rng);
     538                psFree(stats);
     539                return false;
     540            } else {
     541                // Desperate retry
     542                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
     543                stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     544            }
     545        } else {
     546            // Retry with all available pixels
     547            mean = psStatsGetValue(stats, MEAN_STAT);
     548            stdev = psStatsGetValue(stats, STDEV_STAT);
     549        }
     550    } else {
     551        // First attempt
     552        mean = psStatsGetValue(stats, MEAN_STAT);
     553        stdev = psStatsGetValue(stats, STDEV_STAT);
     554    }
     555    psFree(rng);
     556    psFree(stats);
     557    if (!isfinite(mean) || !isfinite(stdev)) {
     558        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     559                "Mean (%f) or stdev (%f) is non-finite.", mean, stdev);
     560        return false;
     561    }
     562
     563    psImage *copy;
     564
     565    switch (image->type.type) {
     566    case PS_TYPE_F32:
     567      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     568      break;
     569    case PS_TYPE_F64:
     570      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     571      break;
     572    default:
     573      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     574      return NULL;
     575      break;
     576    }
     577    // Do scaling
     578    float a = 1.0857362; // 2.5 * log(e);
     579    *bsoften = sqrt(a) * stdev;
     580    //    float m0 = 0; // Can I just arbitrarily set this?
     581   
     582    switch (image->type.type) {
     583    case PS_TYPE_F32:
     584      for (int y = 0; y < numRows; y++) {
     585        for (int x = 0; x < numCols; x++) {
     586/*        if (x == 2331 && y == 2843) { */
     587/*          fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */
     588/*        } */
     589          copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     590        }
     591      }
     592      break;
     593    case PS_TYPE_F64:
     594        for (int y = 0; y < numRows; y++) {
     595          for (int x = 0; x < numCols; x++) {
     596            //      fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset));
     597            copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     598          }
     599        }
     600      break;
     601    default:
     602      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     603      return NULL;
     604      break;
     605    }
     606     
     607    // Do regular scaling on the logarithm image
     608    if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) {
     609      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     610      return false;
     611    }
     612    psFree(copy);
     613    return true;
     614}
     615
     616static bool asinhRange(double *bscale, // Scaling, to return
     617                       double *bzero,  // Zero point, to return
     618                       double *bsoften, // asinh softening parameter, to return
     619                       const psImage *image, // Image to scale
     620                       const psImage *mask,  // Mask image
     621                       psImageMaskType maskVal, // value to mask
     622                       const psFitsOptions *options // FITS options
     623                       )
     624{
     625  psAssert(bscale, "impossible");
     626  psAssert(bzero, "impossible");
     627  psAssert(bsoften, "impossible");
     628  psAssert(image, "impossible");
     629  psAssert(options, "impossible");
     630
     631  psTrace("psLib.fits", 3, "Scaling image by asinh method");
     632  int numCols = image->numCols, numRows = image->numRows; // Size of image
     633 
     634    // Measure the mean and stdev
     635    // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a
     636    // mask.
     637    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     638    psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object
     639    double mean, stdev;                                    // Mean and standard deviation
     640    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
     641        // It could be because the image is entirely masked, in which case we don't want to error
     642        bool good = false;              // Any good pixels?
     643
     644
     645// Find good pixels in an image, by image type
     646#define GOOD_PIXELS_CASE(TYPE) \
     647      case PS_TYPE_##TYPE: \
     648        for (int y = 0; y < image->numRows && !good; y++) { \
     649            for (int x = 0; x < image->numCols && !good; x++) { \
     650                if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \
     651                    continue; \
     652                } \
     653                if (!isfinite(image->data.TYPE[y][x])) { \
     654                    continue; \
     655                } \
     656                good = true; \
     657            } \
     658        } \
     659        break;
     660
     661        switch (image->type.type) {
     662            GOOD_PIXELS_CASE(F32);
     663            GOOD_PIXELS_CASE(F64);
     664          default:
     665            psAbort("Unsupported case: %x", image->type.type);
     666        }
     667
     668        if (!good) {
     669            psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0");
     670            psErrorClear();
     671            *bscale = 1.0;
     672            *bzero = 0.0;
     673            *bsoften = NAN;
     674            psFree(rng);
     675            psFree(stats);
     676            return true;
     677        }
     678
     679        // There are some good pixels in there somewhere; psImageBackground just didn't find them
     680        psLogMsg("psLib.fits", PS_LOG_DETAIL,
     681                 "Couldn't measure background statistics for image quantisation; retrying.");
     682        psErrorClear();
     683        // Retry using all the available pixels
     684        stats->nSubsample = image->numCols * image->numRows + 1;
     685        if (!psImageStats(stats, image, mask, maskVal)) {
     686            psLogMsg("psLib.fits", PS_LOG_DETAIL,
     687                     "Couldn't measure background statistics for image quantisation (attempt 2); retrying.");
     688            psErrorClear();
     689            // Retry with desperate statistic
     690            stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT;
     691            if (!psImageStats(stats, image, mask, maskVal)) {
     692                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image");
     693                psFree(rng);
     694                psFree(stats);
     695                return false;
     696            } else {
     697                // Desperate retry
     698                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
     699                stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     700            }
     701        } else {
     702            // Retry with all available pixels
     703            mean = psStatsGetValue(stats, MEAN_STAT);
     704            stdev = psStatsGetValue(stats, STDEV_STAT);
     705        }
     706    } else {
     707        // First attempt
     708        mean = psStatsGetValue(stats, MEAN_STAT);
     709        stdev = psStatsGetValue(stats, STDEV_STAT);
     710    }
     711    psFree(rng);
     712    psFree(stats);
     713    if (!isfinite(mean) || !isfinite(stdev)) {
     714        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     715                "Mean (%f) or stdev (%f) is non-finite.", mean, stdev);
     716        return false;
     717    }
     718
     719    psImage *copy;
     720
     721    switch (image->type.type) {
     722    case PS_TYPE_F32:
     723      copy = psImageCopy(NULL,image,PS_TYPE_F32);
     724      break;
     725    case PS_TYPE_F64:
     726      copy = psImageCopy(NULL,image,PS_TYPE_F64);
     727      break;
     728    default:
     729      psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type);
     730      return NULL;
     731      break;
     732    }
     733    // Do scaling
     734    float a = 1.0857362; // 2.5 * log(e);
     735    *bsoften = sqrt(a) * stdev;
     736    // float m0 = 0; // Can I just arbitrarily set this?
     737   
     738    switch (image->type.type) {
     739    case PS_TYPE_F32:
     740      for (int y = 0; y < numRows; y++) {
     741        for (int x = 0; x < numCols; x++) {
     742/*        if (x == 2331 && y == 2843) { */
     743/*          fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */
     744/*        } */
     745          copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     746        }
     747      }
     748      break;
     749    case PS_TYPE_F64:
     750        for (int y = 0; y < numRows; y++) {
     751          for (int x = 0; x < numCols; x++) {
     752            //      fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset));
     753            copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     754          }
     755        }
     756      break;
     757    default:
     758      psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type);
     759      return NULL;
     760      break;
     761    }
     762     
     763    // Do regular scaling on the logarithm image
     764    if (!scaleRange(bscale, bzero, copy, options)) {
     765      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     766      return false;
     767    }
     768    psFree(copy);
     769    return true;
     770}
     771
    458772
    459773//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    461775//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    462776
    463 bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image,
     777bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, double *bsoften,
     778                          long *blank, const psImage *image,
    464779                          const psImage *mask, psImageMaskType maskVal, const psFits *fits)
    465780{
    … …  
    521836        }
    522837        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;
     838      case PS_FITS_SCALE_LOG_RANGE:
     839        if (!logscaleRange(bscale,bzero,boffset,image,options)) {
     840          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range");
     841          return false;
     842        }
     843        break;
     844      case PS_FITS_SCALE_LOG_STDEV_POSITIVE:
     845      case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
     846      case PS_FITS_SCALE_LOG_STDEV_BOTH:
     847        if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) {
     848          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     849          return false;
     850        }
     851        break;
     852      case PS_FITS_SCALE_ASINH_RANGE:
     853        if (!asinhRange(bscale,bzero,bsoften,image,mask,maskVal,options)) {
     854          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range");
     855          return false;
     856        }
     857        break;
     858      case PS_FITS_SCALE_ASINH_STDEV_POSITIVE:
     859      case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE:
     860      case PS_FITS_SCALE_ASINH_STDEV_BOTH:
     861        if (!asinhStdev(bscale, bzero,bsoften, image, mask, maskVal, options)) {
     862          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
     863          return false;
     864        }
     865        break;
     866       
     867       
    537868      case PS_FITS_SCALE_MANUAL:
    538869        *bscale = options->bscale;
    … …  
    563894
    564895
    565 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset,
     896psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften,
    566897                            psRandom *rng)
    567898{
    … …  
    626957                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
    627958                value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
     959              }                                                         \
     960              else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)||        \
     961                (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)||        \
     962                  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \
     963                  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
     964                  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) {       \
     965                value = -1.0 * 1.0857362 * (asinh( (IN)->data.INTYPE[y][x] - bsoften)); \
    628966              }                                                         \
    629967                else { \
    … …  
    652990    case PS_TYPE_##INTYPE: { \
    653991        switch (outType) { \
    654             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8); \
    655             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \
    656             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \
    657             SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \
     992          SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);;   \
     993          SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);;  \
     994          SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);;  \
     995          SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);;  \
    658996          default: \
    659997            psAbort("Should be unreachable."); \
    … …  
    6641002    switch (image->type.type) {
    6651003        SCALE_WRITE_IN_CASE(image, F32, out);
    666         SCALE_WRITE_IN_CASE(image, F64, out);
     1004        SCALE_WRITE_IN_CASE(image, F64, out); 
    6671005      default:
    6681006        psAbort("Should be unreachable.");
    … …  
    6791017// the present time, since cfitsio should apply the scaling itself in the process of reading.  However, we may
    6801018// later desire it (e.g., if we ever make our own FITS implementation).
    681 psImage *psFitsScaleFromDisk(const psImage *image, double boffset)
     1019psImage *psFitsScaleFromDisk(const psImage *image, double boffset, double bsoften)
    6821020{
    6831021    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    7181056      for (int y = 0; y < numRows; y++) { \
    7191057          for (int x = 0; x < numCols; x++) { \
    720             out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \
     1058            if (boffset) { \
     1059              out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \
     1060             }                                                                  \
     1061            else if (bsoften) { \
     1062              out->data.OUTTYPE[y][x] = sinh(image->data.INTYPE[y][x] * 2 * bsoften) / (-1.0 * 1.0857362);; \
     1063            } \
    7211064          } \
    7221065      } \
    … …  
    7661109    if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE;
    7671110    if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH;
     1111    if (strcasecmp(string, "ASINH_RANGE") == 0)      return PS_FITS_SCALE_ASINH_RANGE;
     1112    if (strcasecmp(string, "ASINH_MANUAL") == 0)      return PS_FITS_SCALE_ASINH_MANUAL;
     1113    if (strcasecmp(string, "ASINH_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_POSITIVE;
     1114    if (strcasecmp(string, "ASINH_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_NEGATIVE;
     1115    if (strcasecmp(string, "ASINH_STDEV_BOTH") == 0) return PS_FITS_SCALE_ASINH_STDEV_BOTH;
    7681116    if (strcasecmp(string, "MANUAL") == 0)         return PS_FITS_SCALE_MANUAL;
    7691117

  • branches/czw_branch/20110406/psLib/src/fits/psFitsScale.h

    r30636 r31618  
    1111                          double *bzero, ///< Zero point, to return
    1212                          double *boffset, ///< Log offset, to return
     13                          double *bsoften, ///< asinh softening parameter, to return
    1314                          long *blank,  ///< Blank value, to return
    1415                          const psImage *image, ///< Image to scale
    … …  
    2930                            double bzero, ///< Zero point
    3031                            double boffset, ///< Log offset
     32                            double bsoften, ///< asinh softening parameter
    3133                            psRandom *rng ///< Random number generator (for the "fuzz"), or NULL
    3234    );
    3335psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET
    34                              double boffset        ///< Log offset
     36                             double boffset,        ///< Log offset
     37                             double bsoften         ///< asinh softening parameter
    3538                             );
    3639/// Interpret a string as a scaling method

  • branches/czw_branch/20110406/psModules/src/camera/pmFPAfileDefine.c

    r30636 r31618  
    283283              case PS_FITS_SCALE_RANGE:
    284284              case PS_FITS_SCALE_LOG_RANGE:
     285              case PS_FITS_SCALE_ASINH_RANGE:
    285286                // No options required
    286287                break;
    … …  
    289290              case PS_FITS_SCALE_LOG_STDEV_POSITIVE:
    290291              case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
     292              case PS_FITS_SCALE_ASINH_STDEV_POSITIVE:
     293              case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE:
    291294                options->stdevNum = parseOptionFloat(scheme, "STDEV.NUM", source); // Padding to edge
    292295                if (!isfinite(options->stdevNum)) {
    … …  
    299302              case PS_FITS_SCALE_STDEV_BOTH:
    300303              case PS_FITS_SCALE_LOG_STDEV_BOTH:
     304              case PS_FITS_SCALE_ASINH_STDEV_BOTH:
    301305                options->stdevBits = parseOptionInt(scheme, "STDEV.BITS", source, 0); // Bits for stdev
    302306                if (options->stdevBits <= 0) {
    … …  
    312316                options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point
    313317                break;
    314             case PS_FITS_SCALE_LOG_MANUAL:
    315               options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling
    316               options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point
    317               options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset
    318               break;         
    319               default:
     318              case PS_FITS_SCALE_LOG_MANUAL:
     319                options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling
     320                options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point
     321                options->boffset = parseOptionDouble(scheme, "BOFFSET", source); // Log offset
     322              case PS_FITS_SCALE_ASINH_MANUAL:
     323                options->bscale = parseOptionDouble(scheme, "BSCALE", source); // Scaling
     324                options->bzero = parseOptionDouble(scheme, "BZERO", source); // Zero point
     325                options->bsoften = parseOptionDouble(scheme, "BSOFTEN", source); // Softening parameter
     326                break;       
     327              default:
    320328                psAbort("Should never get here.");
    321329            }

  • branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c

    r31607 r31618  
    2222#include <pslib.h>
    2323
     24#include <gsl/gsl_cdf.h>
     25
    2426#include "pmHDU.h"
    2527#include "pmFPA.h"
    … …  
    3537
    3638
    37 # if (0)
     39# if (1)
    3840#define TESTING                         // Enable test output
    39 #define TEST_X 3145                       // x coordinate to examine
    40 #define TEST_Y 2334                       // y coordinate to examine
     41#define TEST_X 5745                       // x coordinate to examine
     42#define TEST_Y 5331                       // y coordinate to examine
    4143#define TEST_RADIUS 0.5                 // Radius to examine
    4244# endif
    … …  
    332334    memset (nGoodBits, 0, 16*sizeof(int));
    333335
    334     // Extract the pixel and mask data
     336
     337    // Extract the pixel and mask data   
    335338    int numGood = 0;                    // Number of good pixels
    336339    for (int i = 0, j = 0; i < inputs->n; i++) {
    … …  
    364367            continue;
    365368        }
    366 
     369       
    367370        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
    368371        psImage *mask = data->readout->mask; // Mask of interest
    … …  
    984987
    985988static void KMMcalculate(const psVector *values,
    986                          float *Punimodal,
     989                         float *Punimodal,int *iter,
    987990                         float *pi1, float *m1, float *s1,
    988991                         float *pi2, float *m2, float *s2) {
    … …  
    10121015  float dL = 0;
    10131016  float oldL = -999;
    1014   int k = 0;
     1017  *iter = 0;
    10151018  logL_bimodal = logL_unimodal;
    10161019  *m1 = mU - 3 * sU;
    … …  
    10241027  float w1,w2;
    10251028
    1026   while (((dL > KMM_TOLERANCE)||(k < 3))&&(k < KMM_MAX_ITERATIONS)) {
    1027     k++;
     1029  float KMM_TOLERANCE = 1e-6;
     1030  int KMM_MAX_ITERATIONS = 500;
     1031  float KMM_SMALL_NUMBER = 1e-5;
     1032  while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
     1033    *iter += 1;
    10281034    dL = fabs(logL_bimodal - oldL);
    10291035    oldL = logL_bimodal;
    … …  
    10551061      *m2 += values->data.F32[i] * P2->data.F32[i];
    10561062
    1057       w1 += P1[i];
    1058       w2 += P2[i];
     1063      w1 += P1->data.F32[i];
     1064      w2 += P2->data.F32[i];
    10591065    }
    10601066    *m1 /= w1;
    … …  
    10951101}
    10961102
    1097 static void KMMrejectUnpopular(const psVector *values, psArray *reject) {
     1103static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
     1104  float KMM_MINIMUM_PVALUE = 0.05;
    10981105  float Punimodal,pi1,m1,s1,pi2,m2,s2;
    1099   KMMcalculate(values,&Punimodal,
     1106  int iter;
     1107  int j;
     1108
     1109  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
     1110  for (j = 0; j < inputs->n; j++) {
     1111    pmStackData *data = inputs->data[j]; // Stack data of interest
     1112    psImage *image = data->readout->image; // Image of interest
     1113    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
     1114    values->data.F32[j] = image->data.F32[yIn][xIn];
     1115  }
     1116 
     1117  KMMcalculate(values,&Punimodal,&iter,
    11001118               &pi1,&m1,&s1,
    11011119               &pi2,&m2,&s2);
     1120  //  fprintf(stderr,
     1121  psTrace("psModules.imcombine",3,
     1122          "KMM Unpopular Test: %d,%d: Puni: %f in %d",x,y,Punimodal,iter);
     1123  //  CHECKPIX(x, y, "
     1124 
    11021125  if (Punimodal < KMM_MINIMUM_PVALUE) {
    11031126    int i;
    1104     float g1,g2;
     1127    float g1,g2,norm;
    11051128    float P1,P2;
    11061129
    … …  
    11121135      P2 = (pi2 * g2) / norm;
    11131136
     1137      CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
     1138               x, y, i, values->data.F32[i],
     1139               P1,m1,s1,pi1,
     1140               P2,m2,s2,pi2,
     1141               (pi1 > pi2)&&(P1 < P2),
     1142               (pi1 < pi2)&&(P1 > P2));
    11141143      if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
    1115         reject_input(reject,i);
     1144        combineMarkReject(inputs,x,y,i);
    11161145      }
    11171146      if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
    1118         reject_input(reject,i);
     1147        combineMarkReject(inputs,x,y,i);
    11191148      }
    11201149    }
    11211150  }
     1151  psFree(values);
    11221152  // else do nothing.
    11231153}
    11241154
    1125 static void KMMrejectBright(const psVector *values, psArray *reject) {
    1126   KMMcalculate(values,&Punimodal,
     1155static void KMMRejectBright(const psArray *inputs, int x, int y) {
     1156  float KMM_MINIMUM_PVALUE = 0.05;
     1157  float Punimodal,pi1,m1,s1,pi2,m2,s2;
     1158  int iter;
     1159  int j;
     1160
     1161  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
     1162  for (j = 0; j < inputs->n; j++) {
     1163    pmStackData *data = inputs->data[j]; // Stack data of interest
     1164    psImage *image = data->readout->image; // Image of interest
     1165    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
     1166    values->data.F32[j] = image->data.F32[yIn][xIn];
     1167  }
     1168 
     1169  KMMcalculate(values,&Punimodal,&iter,
    11271170               &pi1,&m1,&s1,
    11281171               &pi2,&m2,&s2);
    11291172  if (Punimodal < KMM_MINIMUM_PVALUE) {
    11301173    int i;
    1131     float g1,g2;
     1174    float g1,g2,norm;
    11321175    float P1,P2;
    11331176
    … …  
    11401183
    11411184      if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
    1142         reject_input(reject,i);
     1185        combineMarkReject(inputs,x,y,i);
    11431186      }
    11441187      if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
    1145         reject_input(reject,i);
     1188        combineMarkReject(inputs,x,y,i);
    11461189      }
    11471190    }
    … …  
    13101353    }
    13111354
     1355    // Pre-reject inputs using KMM bimodality test.
     1356    if (1)  {
     1357/*       KMMRejectUnpopular(input,x,y); */
     1358      rejection = true;
     1359    }
     1360   
    13121361    // Set up rejection list
    13131362    psArray *pixelMap = NULL;           // Map of pixels to source
    … …  
    13191368    for (int y = minInputRows; y < maxInputRows; y++) {
    13201369        for (int x = minInputCols; x < maxInputCols; x++) {
     1370
    13211371            CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
     1372
     1373            // Pre-reject inputs using KMM bimodality test.
     1374          if (1)  {
     1375            KMMRejectUnpopular(input,x,y);
     1376/*          rejection = true; */
     1377          }
     1378          else {
     1379            KMMRejectBright(input,x,y);
     1380          }
    13221381
    13231382#ifdef TESTING
Note: See TracChangeset for help on using the changeset viewer.