Index: /branches/eam_branches/ipp-20110505/ippScripts/scripts/nightly_science.pl =================================================================== --- /branches/eam_branches/ipp-20110505/ippScripts/scripts/nightly_science.pl (revision 31655) +++ /branches/eam_branches/ipp-20110505/ippScripts/scripts/nightly_science.pl (revision 31656) @@ -60,5 +60,5 @@ 'test_mode' => \$test_mode, 'force_stack_count' => \$force_stack_count, - 'force_diff_count' => \$force_diff_count, + 'force_diff_count' => \$force_diff_count, 'this_target_only=s' => \$this_target_only, 'this_filter_only=s' => \$this_filter_only, @@ -865,47 +865,47 @@ # } -# sub get_lunation_extent { -# my $date = shift; -# my ($year,$month,$day) = split /-/,$date; -# my $dateobs_begin; -# my $dateobs_end; - -# my $dt = DateTime->new(year => $year, month => $month, day => $day, -# hour => 0, minute => 0, second => 0, nanosecond => 0, -# time_zone => 'Pacific/Honolulu'); -# do { -# $dt->subtract(days => 1); -# my $ymd = $dt->ymd; -# my $md_cmd = "moondata $ymd 0 0"; -# my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = -# run ( command => $md_cmd, verbose => $verbose); -# unless ($success) { -# $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); -# &my_die("Unable to perform moondata: $error_code", 0,0, $date, $PS_EXIT_SYS_ERROR); -# } -# my @result = split /\s+/,(join "\n", @$stdout_buf); -# if (abs($result[6]) <= 0.5) { -# $dateobs_end = $ymd; -# } -# } while (!defined($dateobs_end)); - -# do { -# $dt->subtract(days => 1); -# my $ymd = $dt->ymd; -# my $md_cmd = "moondata $ymd 0 0"; -# my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = -# run ( command => $md_cmd, verbose => $verbose); -# unless ($success) { -# $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); -# &my_die("Unable to perform moondata: $error_code", 0,0, $date, $PS_EXIT_SYS_ERROR); -# } -# my @result = split /\s+/,(join "", @$stdout_buf); -# if (abs($result[6]) <= 0.5) { -# $dateobs_begin = $ymd; -# } -# } while (!defined($dateobs_begin)); +sub get_lunation_extent { + my $date = shift; + my ($year,$month,$day) = split /-/,$date; + my $dateobs_begin; + my $dateobs_end; + + my $dt = DateTime->new(year => $year, month => $month, day => $day, + hour => 0, minute => 0, second => 0, nanosecond => 0, + time_zone => 'Pacific/Honolulu'); + do { + $dt->add(days => 1); + my $ymd = $dt->ymd; + my $md_cmd = "moondata $ymd 0 0"; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run ( command => $md_cmd, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform moondata: $error_code", 0,0, $date, $PS_EXIT_SYS_ERROR); + } + my @result = split /\s+/,(join "\n", @$stdout_buf); + if (abs($result[6]) <= 0.6) { + $dateobs_end = $ymd; + } + } while (!defined($dateobs_end)); + + do { + $dt->subtract(days => 1); + my $ymd = $dt->ymd; + my $md_cmd = "moondata $ymd 0 0"; + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run ( command => $md_cmd, verbose => $verbose); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform moondata: $error_code", 0,0, $date, $PS_EXIT_SYS_ERROR); + } + my @result = split /\s+/,(join "", @$stdout_buf); + if (abs($result[6]) <= 0.6) { + $dateobs_begin = $ymd; + } + } while (!defined($dateobs_begin)); -# return($dateobs_begin,$dateobs_end); -# } + return($dateobs_begin,$dateobs_end); +} # sub pre_sweetspot_queue { @@ -1405,6 +1405,12 @@ } } - } - } + if ((defined($science_config{$target}{OFFNIGHT_DIFFS})) && ($science_config{$target}{OFFNIGHT_DIFFS} == 1)) { + foreach my $filter (@filter_list) { + offnight_diff_queue($date,$target,$filter,$pretend); + } + } + } + } + # if (($Npotential == $Nnoexp)&&($metadata_out{nsObservingState} eq 'END_OF_NIGHT')&&($is_processing == 0)) { # $metadata_out{nsDiffState} = 'FINISHED_DIFFS'; @@ -1682,4 +1688,43 @@ } + +sub offnight_diff_queue { + my $date = shift; + my $target = shift; + my $filter = shift; + my $pretend = shift; + my ($label,$workdir,$obs_mode,$object,$comment,$tess_id,$dist_group,$data_group,$reduction) = get_tool_parameters($date,$target); + + my $db = init_gpc_db(); + +# my $Npotential = 0; +# my $Nqueued = 0; + + my ($lunation_start,$lunation_end) = get_lunation_extent($date); + + my $new_data_group = "${data_group}.offnight"; + my $new_dist_group = "${dist_group}.offnight"; + my $cmd = "$difftool -dbname $dbname -definewarpwarp "; + $cmd .= "-input_label $label -template_label $label "; + $cmd .= "-set_workdir $workdir -set_dist_group $new_dist_group -set_data_group $new_data_group "; + $cmd .= " -simple -set_label $label -filter $filter "; + $cmd .= " -dateobs_begin ${lunation_start}T00:00:00 -dateobs_end ${lunation_end}T23:59:59 -distance 1.5 "; + + if (defined($pretend)) { + $cmd .= ' -pretend '; + } + if (defined($debug)) { + print STDERR "ON_diffs wants to run this command: $cmd\n"; + } + + if (($debug == 0)&&(!defined($pretend))) { + my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) = + run ( command => $cmd, verbose => $verbose ); + unless ($success) { + $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR); + &my_die("Unable to perform difftool: $error_code", 0,0,$date, $PS_EXIT_SYS_ERROR); + } + } +} Index: /branches/eam_branches/ipp-20110505/ippTasks/diskbalance.pro =================================================================== --- /branches/eam_branches/ipp-20110505/ippTasks/diskbalance.pro (revision 31655) +++ /branches/eam_branches/ipp-20110505/ippTasks/diskbalance.pro (revision 31656) @@ -93,11 +93,12 @@ $BALANCE_SO_ID_START = int($BALANCE_BIG_SO_ID * rnd(0)) $BALANCE_SO_ID_RANGE = 500000 - +$BALANCE_LIMIT = 500 +$BALANCE_OFFSET = 0 # Select Nebulous objects which should be shuffled task balance.load host local - periods -poll 0.5 - periods -exec 5 + periods -poll 60 + periods -exec 60 periods -timeout 1500 npending 1 @@ -117,5 +118,5 @@ book getword balanceControl control mDestinations -var BALANCE_M - command neb-admin --host $NEB_HOST --db $NEB_DB --user $NEB_USER --pass $NEB_PASS --pendingbalance --limit 500 --so_id_start $BALANCE_SO_ID_START --so_id_range $BALANCE_SO_ID_RANGE --balancesources $BALANCE_N --balancedestination $BALANCE_M + command neb-admin --host $NEB_HOST --db $NEB_DB --user $NEB_USER --pass $NEB_PASS --pendingbalance --limit $BALANCE_LIMIT --offset $BALANCE_OFFSET --so_id_start $BALANCE_SO_ID_START --so_id_range $BALANCE_SO_ID_RANGE --balancesources $BALANCE_N --balancedestination $BALANCE_M end @@ -125,5 +126,5 @@ # advance the so_id counter $BALANCE_SO_ID_START = $BALANCE_SO_ID_START + $BALANCE_SO_ID_RANGE - + $BALANCE_OFFSET = 0 # convert 'stdout' to book format ipptool2book stdout balancePending -key key -uniq -setword pantaskState INIT @@ -141,5 +142,5 @@ # convert 'stdout' to book format ipptool2book stdout balancePending -key key -uniq -setword pantaskState INIT - + $BALANCE_OFFSET = $BALANCE_OFFSET + $BALANCE_LIMIT if ($VERBOSE > 2) book listbook balancePending Index: /branches/eam_branches/ipp-20110505/ippTools/src/stacktool.c =================================================================== --- /branches/eam_branches/ipp-20110505/ippTools/src/stacktool.c (revision 31655) +++ /branches/eam_branches/ipp-20110505/ippTools/src/stacktool.c (revision 31656) @@ -243,6 +243,6 @@ PXOPT_COPY_STR(config->args, where, "-select_exp_type", "rawExp.exp_type", "=="); PXOPT_COPY_F32(config->args, where, "-select_good_frac_min", "warpSkyfile.good_frac", ">="); - PXOPT_COPY_STR(config->args, where, "-select_skycell_id", "warpSkyfile.skycell_id", "LIKE"); - PXOPT_COPY_STR(config->args, where, "-select_data_group", "warpRun.data_group", "=="); + pxAddLabelSearchArgs(config, where, "-select_skycell_id", "warpSkyfile.skycell_id", "LIKE"); + pxAddLabelSearchArgs(config, where, "-select_data_group", "warpRun.data_group", "LIKE"); pxAddLabelSearchArgs (config, where, "-select_label", "warpRun.label", "LIKE"); // define using warp label pxAddLabelSearchArgs (config, where, "-warp_id", "warpRun.warp_id", "=="); Index: /branches/eam_branches/ipp-20110505/ippTools/src/stacktoolConfig.c =================================================================== --- /branches/eam_branches/ipp-20110505/ippTools/src/stacktoolConfig.c (revision 31655) +++ /branches/eam_branches/ipp-20110505/ippTools/src/stacktoolConfig.c (revision 31656) @@ -55,9 +55,9 @@ psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-set_dvodb", 0, "define dvodb", NULL); psMetadataAddTime(definebyqueryArgs, PS_LIST_TAIL, "-set_registered", 0, "time detrend run was registered", now); - psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_skycell_id", 0, "search for skycell_id", NULL); + psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_skycell_id", PS_META_DUPLICATE_OK, "search for skycell_id (LIKE comparison, multiple ok)", NULL); psMetadataAddF32(definebyqueryArgs, PS_LIST_TAIL, "-select_good_frac_min", 0, "define min good_frac", 0.0); psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_exp_type", 0, "search for exp_type", "object"); psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_label", PS_META_DUPLICATE_OK, "search by warpRun label (LIKE comparison, multiple OK)", NULL); - psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_data_group", 0 , "search by warpRun data_group", NULL); + psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_data_group", PS_META_DUPLICATE_OK , "search by warpRun data_group (LIKE comparison, multiple ok)", NULL); psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_comment", 0, "search for comment (LIKE)", NULL); psMetadataAddStr(definebyqueryArgs, PS_LIST_TAIL, "-select_inst", 0, "search for camera", NULL); Index: /branches/eam_branches/ipp-20110505/psLib/src/fits/psFits.h =================================================================== --- /branches/eam_branches/ipp-20110505/psLib/src/fits/psFits.h (revision 31655) +++ /branches/eam_branches/ipp-20110505/psLib/src/fits/psFits.h (revision 31656) @@ -55,5 +55,10 @@ PS_FITS_SCALE_LOG_STDEV_NEGATIVE, ///< Take logarithm, Auto-scale to sample stdev, place mean at upper limit PS_FITS_SCALE_LOG_STDEV_BOTH, ///< Take logarithm, Auto-scale to sample stdev, place mean at middle - PS_FITS_SCALE_LOG_MANUAL ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) + PS_FITS_SCALE_LOG_MANUAL, ///< Manual scaling (use specified BSCALE, BZERO, and BOFFSET) + PS_FITS_SCALE_ASINH_RANGE, ///< Do asinh scaling, auto-scale to preserve dynamic range + PS_FITS_SCALE_ASINH_STDEV_POSITIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at lower limit + PS_FITS_SCALE_ASINH_STDEV_NEGATIVE, ///< Do asinh scaling, auto-scale to sample stdev, place mean at upper limit + PS_FITS_SCALE_ASINH_STDEV_BOTH, ///< Do asinh scaling, auto-scale to sample stdev, place mean at middle + PS_FITS_SCALE_ASINH_MANUAL ///< Manual scaling after doing asinh scaling. } psFitsScaling; @@ -72,4 +77,5 @@ double bscale, bzero; ///< Manually specified BSCALE and BZERO (for SCALE_MANUAL) double boffset; ///< Manually specified BOFFSET (for SCALE_MANUAL) + double bsoften; ///< Manuallly specified BSOFTEN (for SCALE_MANUAL) double mean, stdev; ///< Mean and standard deviation of image int stdevBits; ///< Number of bits to sample a standard deviation (for SCALE_STDEV_*) Index: /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsImage.c =================================================================== --- /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsImage.c (revision 31655) +++ /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsImage.c (revision 31656) @@ -55,4 +55,5 @@ int psDatatype; // psLib data type bool is_logscaled; // is this image log scaled using BOFFSET? + bool is_asinh; // is this image asinh scaled using BSOFTEN? } p_psFitsReadInfo; @@ -129,5 +130,5 @@ // Check scale and zero - double bscale = 0.0, bzero = 0.0, boffset = NAN; // Scale and zero point + double bscale = 0.0, bzero = 0.0, boffset = NAN, bsoften = NAN; // Scale and zero point if (fits_read_key_dbl(fits->fd, "BSCALE", &bscale, NULL, &status) && status != KEY_NO_EXIST) { psFitsError(status, true, "Unable to read header."); @@ -153,4 +154,20 @@ info->is_logscaled = false; } + status = 0; + if (fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status) && status != KEY_NO_EXIST) { + psFitsError(status, true, "Unable to read header."); + goto bad; + } + if (status == KEY_NO_EXIST) { + info->is_asinh = false; + } + else if (isfinite(bsoften)) { + info->is_asinh = true; + info->is_logscaled = false; + } + else { + info->is_asinh = false; + } + status = 0; @@ -262,4 +279,5 @@ double *bzero, // Zero point applied double *boffset, // Log offset applied + double *bsoften, // asinh offset applied long *blank, // Blank value (integer data) psFitsFloat *floatType, // Type of custom floating-point @@ -275,4 +293,5 @@ psAssert(bzero, "impossible"); psAssert(boffset, "impossible"); + psAssert(bsoften, "impossible"); psAssert(floatType, "impossible"); psAssert(fits, "impossible"); @@ -322,5 +341,5 @@ if (newScaleZero) { // Choose an appropriate BSCALE and BZERO - if (!psFitsScaleDetermine(bscale, bzero, boffset, blank, image, mask, maskVal, fits)) { + if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) { // We can't have the write dying for this reason --- try to save it somehow! psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise."); @@ -348,5 +367,5 @@ } - return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, rng); + return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng); } @@ -480,6 +499,6 @@ outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType); } - // Need to apply BOFFSET if info->is_logscaled is true + if (info->is_logscaled) { double boffset; @@ -487,5 +506,15 @@ status = 0; fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); - psImage *newImage = psFitsScaleFromDisk(outImage,boffset); + psImage *newImage = psFitsScaleFromDisk(outImage,boffset,0.0); + psFree (outImage); + outImage = newImage; + } + else if (info->is_asinh) { + double bsoften,boffset; + int status; + status = 0; + fits_read_key_dbl(fits->fd, "BSOFTEN", &bsoften, NULL, &status); + fits_read_key_dbl(fits->fd, "BOFFSET", &boffset, NULL, &status); + psImage *newImage = psFitsScaleFromDisk(outImage,boffset,bsoften); psFree (outImage); outImage = newImage; @@ -602,7 +631,8 @@ double bscale = NAN, bzero = NAN; // Scale and zero point to put in header double boffset = NAN; // Log offset to put into header. + double bsoften = NAN; // Asinh softening parameter to put into header long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, image, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image, mask, maskVal, NULL, true); // Image to write out if (!diskImage) { @@ -696,9 +726,24 @@ if (options&&(!((options->scaling == PS_FITS_SCALE_LOG_RANGE)|| (options->scaling == PS_FITS_SCALE_LOG_MANUAL)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| - (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { + (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)|| + (options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { if (psMetadataLookup(header,"BOFFSET")) { psMetadataRemoveKey(header,"BOFFSET"); + } + } + // Remove any BSOFTENvalues that exist in the header if we are not using that scaling anymore + if (options&&(!((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { + if (psMetadataLookup(header,"BSOFTEN")) { + psMetadataRemoveKey(header,"BSOFTEN"); } } @@ -727,6 +772,19 @@ (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)))) { fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, - "Scaling: TRUE = BZERO + BSCALE * 10**(DISK) + BOFFSET)", &status); + "Scaling: ORIGINAL = 10**(TRUE) + BOFFSET)", &status); } + if (options&&(((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)))) { + fits_write_key_dbl(fits->fd, "BOFFSET", boffset, 12, + "Scaling: NORM = DISK + BOFFSET", &status); + fits_write_key_dbl(fits->fd, "BSOFTEN", bsoften, 12, + "Scaling: TRUE = 2 * BSOFTEN * sinh(NORM/a)",&status); + + + } + if (psFitsError(status, true, "Could not write BSCALE/BZERO headers to file.")) { success = false; @@ -831,7 +889,8 @@ double bscale = NAN, bzero = NAN; // Scale and zero point to put in header double boffset = NAN; // Log offset to put in header + double bsoften = NAN; // Asinh softening parameter to put in header long blank = 0; // Blank (undefined) value for image psFitsFloat floatType; // Custom floating-point convention type - psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &blank, &floatType, fits, input, + psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input, mask, maskVal, NULL, false); // Image to write out if (!diskImage) { Index: /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.c =================================================================== --- /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.c (revision 31655) +++ /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.c (revision 31656) @@ -214,4 +214,5 @@ case PS_FITS_SCALE_STDEV_POSITIVE: case PS_FITS_SCALE_LOG_STDEV_POSITIVE: + case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: // Put (mean - N sigma) at the lowest possible value: predominantly positive images imageVal = mean - options->stdevNum * stdev; @@ -220,4 +221,5 @@ case PS_FITS_SCALE_STDEV_NEGATIVE: case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: + case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: // Put (mean + N sigma) at the highest possible value: predominantly negative images imageVal = mean + options->stdevNum * stdev; @@ -226,4 +228,5 @@ case PS_FITS_SCALE_STDEV_BOTH: case PS_FITS_SCALE_LOG_STDEV_BOTH: + case PS_FITS_SCALE_ASINH_STDEV_BOTH: // Put mean right in the middle: images with an equal abundance of positive and negative values imageVal = mean; @@ -241,4 +244,6 @@ } + + static bool logscaleStdev(double *bscale, // Scaling, to return @@ -456,4 +461,347 @@ } +static bool asinhStdev(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // asinh flux zeropoint, to return + double *bsoften, // asinh softening parameter, to return + const psImage *image, // Image to scale + const psImage *mask, // Mask image + psImageMaskType maskVal, // value to mask + const psFitsOptions *options // FITS options + ) +{ + psAssert(bscale, "impossible"); + psAssert(bzero, "impossible"); + psAssert(boffset, "impossible"); + psAssert(bsoften, "impossible"); + psAssert(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by asinh method"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Measure the mean and stdev + // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a + // mask. + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + double mean, stdev; // Mean and standard deviation + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + // It could be because the image is entirely masked, in which case we don't want to error + bool good = false; // Any good pixels? + + +// Find good pixels in an image, by image type +#define GOOD_PIXELS_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int y = 0; y < image->numRows && !good; y++) { \ + for (int x = 0; x < image->numCols && !good; x++) { \ + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ + continue; \ + } \ + if (!isfinite(image->data.TYPE[y][x])) { \ + continue; \ + } \ + good = true; \ + } \ + } \ + break; + + switch (image->type.type) { + GOOD_PIXELS_CASE(F32); + GOOD_PIXELS_CASE(F64); + default: + psAbort("Unsupported case: %x", image->type.type); + } + + if (!good) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); + psErrorClear(); + *bscale = 1.0; + *bzero = 0.0; + *bsoften = NAN; + psFree(rng); + psFree(stats); + return true; + } + + // There are some good pixels in there somewhere; psImageBackground just didn't find them + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation; retrying."); + psErrorClear(); + // Retry using all the available pixels + stats->nSubsample = image->numCols * image->numRows + 1; + if (!psImageStats(stats, image, mask, maskVal)) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); + psErrorClear(); + // Retry with desperate statistic + stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; + if (!psImageStats(stats, image, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); + psFree(rng); + psFree(stats); + return false; + } else { + // Desperate retry + mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); + stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); + } + } else { + // Retry with all available pixels + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + } else { + // First attempt + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + psFree(rng); + psFree(stats); + if (!isfinite(mean) || !isfinite(stdev)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); + return false; + } + + psImage *copy; + + switch (image->type.type) { + case PS_TYPE_F32: + copy = psImageCopy(NULL,image,PS_TYPE_F32); + break; + case PS_TYPE_F64: + copy = psImageCopy(NULL,image,PS_TYPE_F64); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); + return NULL; + break; + } + // Do scaling + float a = 1.0857362; // 2.5 * log(e); + *bsoften = sqrt(a) * stdev; + *boffset = mean; + // float m0 = 0; // Can I just arbitrarily set this? + + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { +/* if (x == 266 && y == 4584) { */ +/* 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), */ +/* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ +/* } */ + if (isfinite(image->data.F32[y][x])) { + copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F32[y][x] = image->data.F32[y][x]; + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); + if (isfinite(image->data.F64[y][x])) { + copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F64[y][x] = image->data.F64[y][x]; + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + + // Do regular scaling on the logarithm image + if (!scaleStdev(bscale, bzero, copy, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + psFree(copy); + return true; +} + +static bool asinhRange(double *bscale, // Scaling, to return + double *bzero, // Zero point, to return + double *boffset, // asinh flux zeropoint, to return + double *bsoften, // asinh softening parameter, to return + const psImage *image, // Image to scale + const psImage *mask, // Mask image + psImageMaskType maskVal, // value to mask + const psFitsOptions *options // FITS options + ) +{ + psAssert(bscale, "impossible"); + psAssert(bzero, "impossible"); + psAssert(bsoften, "impossible"); + psAssert(image, "impossible"); + psAssert(options, "impossible"); + + psTrace("psLib.fits", 3, "Scaling image by asinh method"); + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Measure the mean and stdev + // psImageBackground automatically excludes pixels that are non-finite, so we don't need to bother about a + // mask. + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object + double mean, stdev; // Mean and standard deviation + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + // It could be because the image is entirely masked, in which case we don't want to error + bool good = false; // Any good pixels? + + +// Find good pixels in an image, by image type +#define GOOD_PIXELS_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + for (int y = 0; y < image->numRows && !good; y++) { \ + for (int x = 0; x < image->numCols && !good; x++) { \ + if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \ + continue; \ + } \ + if (!isfinite(image->data.TYPE[y][x])) { \ + continue; \ + } \ + good = true; \ + } \ + } \ + break; + + switch (image->type.type) { + GOOD_PIXELS_CASE(F32); + GOOD_PIXELS_CASE(F64); + default: + psAbort("Unsupported case: %x", image->type.type); + } + + if (!good) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0"); + psErrorClear(); + *bscale = 1.0; + *bzero = 0.0; + *bsoften = NAN; + psFree(rng); + psFree(stats); + return true; + } + + // There are some good pixels in there somewhere; psImageBackground just didn't find them + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation; retrying."); + psErrorClear(); + // Retry using all the available pixels + stats->nSubsample = image->numCols * image->numRows + 1; + if (!psImageStats(stats, image, mask, maskVal)) { + psLogMsg("psLib.fits", PS_LOG_DETAIL, + "Couldn't measure background statistics for image quantisation (attempt 2); retrying."); + psErrorClear(); + // Retry with desperate statistic + stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT; + if (!psImageStats(stats, image, mask, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image"); + psFree(rng); + psFree(stats); + return false; + } else { + // Desperate retry + mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT); + stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT); + } + } else { + // Retry with all available pixels + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + } else { + // First attempt + mean = psStatsGetValue(stats, MEAN_STAT); + stdev = psStatsGetValue(stats, STDEV_STAT); + } + psFree(rng); + psFree(stats); + if (!isfinite(mean) || !isfinite(stdev)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Mean (%f) or stdev (%f) is non-finite.", mean, stdev); + return false; + } + + psImage *copy; + + switch (image->type.type) { + case PS_TYPE_F32: + copy = psImageCopy(NULL,image,PS_TYPE_F32); + break; + case PS_TYPE_F64: + copy = psImageCopy(NULL,image,PS_TYPE_F64); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Target type is not a float: %d", image->type.type); + return NULL; + break; + } + // Do scaling + float a = 1.0857362; // 2.5 * log(e); + *bsoften = sqrt(a) * stdev; + *boffset = mean; + // float m0 = 0; // Can I just arbitrarily set this? + + switch (image->type.type) { + case PS_TYPE_F32: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { +/* if (x == 2331 && y == 2843) { */ +/* fprintf(stderr,"psFS32: %d %d %g %g %g\n",x,y,offset,image->data.F32[y][x],log10(image->data.F32[y][x] - offset)); */ +/* } */ +/* if (x == 266 && y == 4584) { */ +/* 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), */ +/* a * asinh((image->data.F32[y][x] - *boffset) / (2 * *bsoften))); */ +/* } */ + + if (isfinite(image->data.F32[y][x])) { + copy->data.F32[y][x] = a * asinh( (image->data.F32[y][x] - *boffset) / (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F32[y][x] = image->data.F32[y][x]; + } + } + } + break; + case PS_TYPE_F64: + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // fprintf(stderr,"psFS64: %d %d %g %g %g\n",x,y,offset,image->data.F64[y][x],log10(image->data.F64[y][x] - offset)); + if (isfinite(image->data.F64[y][x])) { + copy->data.F64[y][x] = a * asinh( (image->data.F64[y][x] - *boffset)/ (2 * *bsoften));// - 2.5 * log10(b) + m0; + } + else { + copy->data.F64[y][x] = image->data.F32[y][x]; + } + } + } + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Target type is not a float: %d",image->type.type); + return NULL; + break; + } + + // Do regular scaling on the asinh image + if (!scaleRange(bscale, bzero, copy, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + psFree(copy); + return true; +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -461,5 +809,6 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, long *blank, const psImage *image, +bool psFitsScaleDetermine(double *bscale, double *bzero, double *boffset, double *bsoften, + long *blank, const psImage *image, const psImage *mask, psImageMaskType maskVal, const psFits *fits) { @@ -521,18 +870,34 @@ } break; - case PS_FITS_SCALE_LOG_RANGE: - if (!logscaleRange(bscale,bzero,boffset,image,options)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); - return false; - } - break; - case PS_FITS_SCALE_LOG_STDEV_POSITIVE: - case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: - case PS_FITS_SCALE_LOG_STDEV_BOTH: - if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); - return false; - } - break; + case PS_FITS_SCALE_LOG_RANGE: + if (!logscaleRange(bscale,bzero,boffset,image,options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_LOG_STDEV_POSITIVE: + case PS_FITS_SCALE_LOG_STDEV_NEGATIVE: + case PS_FITS_SCALE_LOG_STDEV_BOTH: + if (!logscaleStdev(bscale, bzero,boffset, image, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + case PS_FITS_SCALE_ASINH_RANGE: + if (!asinhRange(bscale,bzero,boffset,bsoften,image,mask,maskVal,options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from range"); + return false; + } + break; + case PS_FITS_SCALE_ASINH_STDEV_POSITIVE: + case PS_FITS_SCALE_ASINH_STDEV_NEGATIVE: + case PS_FITS_SCALE_ASINH_STDEV_BOTH: + if (!asinhStdev(bscale, bzero,boffset,bsoften, image, mask, maskVal, options)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set BSCALE and BZERO from stdev"); + return false; + } + break; + + case PS_FITS_SCALE_MANUAL: *bscale = options->bscale; @@ -557,11 +922,11 @@ - psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BLANK = %ld\n", - *bscale, *bzero, *boffset, *blank); + psTrace("psLib.fits", 3, "BSCALE = %.10lf, BZERO = %.10lf, BOFFSET = %.10lf, BSOFTEN = %.10lf, BLANK = %ld\n", + *bscale, *bzero, *boffset, *bsoften, *blank); return true; } -psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, +psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, psRandom *rng) { @@ -625,5 +990,26 @@ (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \ (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \ - value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + if (isfinite((IN)->data.INTYPE[y][x])) { \ + value = log10( (IN)->data.INTYPE[y][x] - boffset ); \ + } \ + else { \ + value = (IN)->data.INTYPE[y][x]; \ + } \ + } \ + else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \ + (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \ + if (isfinite((IN)->data.INTYPE[y][x])) { \ + value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \ + } \ + else { \ + value = (IN)->data.INTYPE[y][x]; \ + } \ + if ((x == 1000)&&(y == 1000)) { \ + fprintf(stderr,"ASINH: %f %f %f", \ + (IN)->data.INTYPE[y][x],value,bsoften); \ + } \ } \ else { \ @@ -652,8 +1038,8 @@ case PS_TYPE_##INTYPE: { \ switch (outType) { \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32); \ - SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64); \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);; \ + SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);; \ default: \ psAbort("Should be unreachable."); \ @@ -664,5 +1050,5 @@ switch (image->type.type) { SCALE_WRITE_IN_CASE(image, F32, out); - SCALE_WRITE_IN_CASE(image, F64, out); + SCALE_WRITE_IN_CASE(image, F64, out); default: psAbort("Should be unreachable."); @@ -679,5 +1065,5 @@ // the present time, since cfitsio should apply the scaling itself in the process of reading. However, we may // later desire it (e.g., if we ever make our own FITS implementation). -psImage *psFitsScaleFromDisk(const psImage *image, double boffset) +psImage *psFitsScaleFromDisk(const psImage *image, double boffset, double bsoften) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -718,5 +1104,10 @@ for (int y = 0; y < numRows; y++) { \ for (int x = 0; x < numCols; x++) { \ - out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ + if (bsoften) { \ + out->data.OUTTYPE[y][x] = 2 * bsoften * sinh(image->data.INTYPE[y][x] / (1.0857362)) + boffset; \ + } \ + else if (boffset) { \ + out->data.OUTTYPE[y][x] = pow(10,image->data.INTYPE[y][x]) + boffset;; \ + } \ } \ } \ @@ -766,4 +1157,9 @@ if (strcasecmp(string, "LOG_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_LOG_STDEV_NEGATIVE; if (strcasecmp(string, "LOG_STDEV_BOTH") == 0) return PS_FITS_SCALE_LOG_STDEV_BOTH; + if (strcasecmp(string, "ASINH_RANGE") == 0) return PS_FITS_SCALE_ASINH_RANGE; + if (strcasecmp(string, "ASINH_MANUAL") == 0) return PS_FITS_SCALE_ASINH_MANUAL; + if (strcasecmp(string, "ASINH_STDEV_POSITIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_POSITIVE; + if (strcasecmp(string, "ASINH_STDEV_NEGATIVE") == 0) return PS_FITS_SCALE_ASINH_STDEV_NEGATIVE; + if (strcasecmp(string, "ASINH_STDEV_BOTH") == 0) return PS_FITS_SCALE_ASINH_STDEV_BOTH; if (strcasecmp(string, "MANUAL") == 0) return PS_FITS_SCALE_MANUAL; Index: /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.h =================================================================== --- /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.h (revision 31655) +++ /branches/eam_branches/ipp-20110505/psLib/src/fits/psFitsScale.h (revision 31656) @@ -11,4 +11,5 @@ double *bzero, ///< Zero point, to return double *boffset, ///< Log offset, to return + double *bsoften, ///< asinh softening parameter, to return long *blank, ///< Blank value, to return const psImage *image, ///< Image to scale @@ -29,8 +30,10 @@ double bzero, ///< Zero point double boffset, ///< Log offset + double bsoften, ///< asinh softening parameter psRandom *rng ///< Random number generator (for the "fuzz"), or NULL ); psImage *psFitsScaleFromDisk(const psImage *image, ///< Image to to unapply BOFFSET - double boffset ///< Log offset + double boffset, ///< Log offset + double bsoften ///< asinh softening parameter ); /// Interpret a string as a scaling method