IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 13, 2011, 3:52:49 PM (15 years ago)
Author:
watersc1
Message:

asinh scaling. Both writing and reading seem to function as expected, so I'm calling it done.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r31618 r31621  
    214214      case PS_FITS_SCALE_STDEV_POSITIVE:
    215215      case PS_FITS_SCALE_LOG_STDEV_POSITIVE:
     216      case PS_FITS_SCALE_ASINH_STDEV_POSITIVE:
    216217        // Put (mean - N sigma) at the lowest possible value: predominantly positive images
    217218        imageVal = mean - options->stdevNum * stdev;
    … …  
    220221      case PS_FITS_SCALE_STDEV_NEGATIVE:
    221222      case PS_FITS_SCALE_LOG_STDEV_NEGATIVE:
     223      case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE:
    222224        // Put (mean + N sigma) at the highest possible value: predominantly negative images
    223225        imageVal = mean + options->stdevNum * stdev;
    … …  
    226228      case PS_FITS_SCALE_STDEV_BOTH:
    227229      case PS_FITS_SCALE_LOG_STDEV_BOTH:
     230      case PS_FITS_SCALE_ASINH_STDEV_BOTH:
    228231        // Put mean right in the middle: images with an equal abundance of positive and negative values
    229232        imageVal = mean;
    … …  
    460463static bool asinhStdev(double *bscale, // Scaling, to return
    461464                       double *bzero,  // Zero point, to return
     465                       double *boffset, // asinh flux zeropoint, to return
    462466                       double *bsoften, // asinh softening parameter, to return
    463467                       const psImage *image, // Image to scale
    … …  
    469473  psAssert(bscale, "impossible");
    470474  psAssert(bzero, "impossible");
     475  psAssert(boffset, "impossible");
    471476  psAssert(bsoften, "impossible");
    472477  psAssert(image, "impossible");
    … …  
    578583    float a = 1.0857362; // 2.5 * log(e);
    579584    *bsoften = sqrt(a) * stdev;
     585    *boffset = mean;
    580586    //    float m0 = 0; // Can I just arbitrarily set this?
    581587   
    … …  
    584590      for (int y = 0; y < numRows; y++) {
    585591        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)); */
     592/*        if (x == 266 && y == 4584) { */
     593/*          fprintf(stderr,"psFS32: %d %d %g %g %g %g %g\n",x,y,*boffset,*bsoften,image->data.F32[y][x],log10(image->data.F32[y][x] - *boffset), */
     594/*                  a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */
    588595/*        } */
    589           copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     596          if (isfinite(image->data.F32[y][x])) {
     597            copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     598          }
     599          else {
     600            copy->data.F32[y][x] = image->data.F32[y][x];
     601          }
    590602        }
    591603      }
    … …  
    595607          for (int x = 0; x < numCols; x++) {
    596608            //      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;
     609            if (isfinite(image->data.F64[y][x])) {
     610              copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     611            }
     612            else {
     613              copy->data.F64[y][x] = image->data.F64[y][x];
     614            }     
    598615          }
    599616        }
    … …  
    616633static bool asinhRange(double *bscale, // Scaling, to return
    617634                       double *bzero,  // Zero point, to return
     635                       double *boffset, // asinh flux zeropoint, to return
    618636                       double *bsoften, // asinh softening parameter, to return
    619637                       const psImage *image, // Image to scale
    … …  
    734752    float a = 1.0857362; // 2.5 * log(e);
    735753    *bsoften = sqrt(a) * stdev;
     754    *boffset = mean;
    736755    // float m0 = 0; // Can I just arbitrarily set this?
    737756   
    … …  
    743762/*          fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */
    744763/*        } */
    745           copy->data.F32[y][x] = -1.0 * a * asinh( image->data.F32[y][x] / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     764/*        if (x == 266 && y == 4584) { */
     765/*          fprintf(stderr,"psFS32: %d %d %g %g %g %g %g\n",x,y,*boffset,*bsoften,image->data.F32[y][x],log10(image->data.F32[y][x] - *boffset), */
     766/*                  a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */
     767/*        } */
     768
     769          if (isfinite(image->data.F32[y][x])) {
     770            copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0;
     771          }
     772          else {
     773            copy->data.F32[y][x] = image->data.F32[y][x];
     774          }
    746775        }
    747776      }
    … …  
    751780          for (int x = 0; x < numCols; x++) {
    752781            //      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;
     782            if (isfinite(image->data.F64[y][x])) {
     783              copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset)/ (2 * *bsoften));// - 2.5 * log10(b) + m0;
     784            }
     785            else {
     786              copy->data.F64[y][x] = image->data.F32[y][x];
     787            }
    754788          }
    755789        }
    … …  
    761795    }
    762796     
    763     // Do regular scaling on the logarithm image
     797    // Do regular scaling on the asinh image
    764798    if (!scaleRange(bscale, bzero, copy, options)) {
    765799      psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
    … …  
    851885        break;
    852886      case PS_FITS_SCALE_ASINH_RANGE:
    853         if (!asinhRange(bscale,bzero,bsoften,image,mask,maskVal,options)) {
     887        if (!asinhRange(bscale,bzero,boffset,bsoften,image,mask,maskVal,options)) {
    854888          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range");
    855889          return false;
    … …  
    859893      case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE:
    860894      case PS_FITS_SCALE_ASINH_STDEV_BOTH:
    861         if (!asinhStdev(bscale, bzero,bsoften, image, mask, maskVal, options)) {
     895        if (!asinhStdev(bscale, bzero,boffset,bsoften, image, mask, maskVal, options)) {
    862896          psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev");
    863897          return false;
    … …  
    888922
    889923
    890     psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n",
    891             *bscale, *bzero, *boffset, *blank);
     924    psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BSOFTEN = %.10lf, BLANK = %ld\n",
     925            *bscale, *bzero, *boffset, *bsoften, *blank);
    892926    return true;
    893927}
    … …  
    956990                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
    957991                  (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
    958                 value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
     992                if (isfinite((IN)->data.INTYPE[y][x])) {                \
     993                  value = log10( (IN)->data.INTYPE[y][x] - boffset );   \
     994                }                                                       \
     995                else {                                                  \
     996                  value = (IN)->data.INTYPE[y][x];                      \
     997                }                                                       \
    959998              }                                                         \
    960999              else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)||        \
    … …  
    9631002                  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
    9641003                  (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) {       \
    965                 value = -1.0 * 1.0857362 * (asinh( (IN)->data.INTYPE[y][x] - bsoften)); \
     1004                if (isfinite((IN)->data.INTYPE[y][x])) {                \
     1005                  value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \
     1006                }                                                       \
     1007                else {                                                  \
     1008                  value = (IN)->data.INTYPE[y][x];                      \
     1009                }                                                       \
     1010                if ((x == 1000)&&(y == 1000)) {                         \
     1011                  fprintf(stderr,"ASINH: %f %f %f",                     \
     1012                          (IN)->data.INTYPE[y][x],value,bsoften);       \
     1013                }                                                       \
    9661014              }                                                         \
    9671015                else { \
    … …  
    10561104      for (int y = 0; y < numRows; y++) { \
    10571105          for (int x = 0; x < numCols; x++) { \
    1058             if (boffset) { \
     1106            if (bsoften) {                                              \
     1107              out->data.OUTTYPE[y][x] = 2 * bsoften * sinh(image->data.INTYPE[y][x] / (1.0857362)) + boffset; \
     1108            }                                                           \
     1109            else if (boffset) {                                 \
    10591110              out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \
    10601111             }                                                                  \
    1061             else if (bsoften) { \
    1062               out->data.OUTTYPE[y][x] = sinh(image->data.INTYPE[y][x] * 2 * bsoften) / (-1.0 * 1.0857362);; \
    1063             } \
    10641112          } \
    10651113      } \
Note: See TracChangeset for help on using the changeset viewer.