Index: branches/eam_branches/20091201/ppStack/src/ppStackMatch.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackMatch.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackMatch.c (revision 26674) @@ -317,11 +317,11 @@ float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor" - bool scale = psMetadataLookupBool(NULL, recipe, "SCALE"); // Scale kernel parameters? - float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling - float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling - float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling + bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", scaleRef, scaleMin, scaleMax); return false; Index: branches/eam_branches/20091201/ppStack/src/ppStackOptions.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackOptions.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackOptions.c (revision 26674) @@ -24,4 +24,5 @@ psFree(options->sourceLists); psFree(options->norm); + psFree(options->sources); psFree(options->cells); psFree(options->kernels); @@ -44,4 +45,5 @@ options->convolve = true; options->matchZPs = true; + options->photometry = false; options->stats = NULL; options->statsFile = NULL; @@ -60,4 +62,5 @@ options->sourceLists = NULL; options->norm = NULL; + options->sources = NULL; options->cells = NULL; options->kernels = NULL; Index: branches/eam_branches/20091201/ppStack/src/ppStackOptions.h =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackOptions.h (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackOptions.h (revision 26674) @@ -10,4 +10,5 @@ bool convolve; // Convolve images? bool matchZPs; // Adjust relative fluxes based on transparency analysis? + bool photometry; // Perform photometry? psMetadata *stats; // Statistics for output FILE *statsFile; // File to which to write statistics @@ -24,4 +25,5 @@ psArray *sourceLists; // Individual lists of sources for matching psVector *norm; // Normalisation for each image + psArray *sources; // Matched sources // Convolve psArray *cells; // Cells for convolved images --- a handle for reading again Index: branches/eam_branches/20091201/ppStack/src/ppStackPhotometry.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackPhotometry.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackPhotometry.c (revision 26674) @@ -19,6 +19,5 @@ psAssert(recipe, "We've thrown an error on this before."); - bool mdok; // Status of MD lookup - if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { + if (!options->photometry) { // Nothing to do return true; @@ -60,6 +59,5 @@ "Bits to use for mark", markValue); - pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT"); - psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES"); + psArray *inSources = options->sources; if (!inSources) { psError(PS_ERR_UNKNOWN, false, "Unable to find input sources"); Index: branches/eam_branches/20091201/ppStack/src/ppStackReadout.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackReadout.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackReadout.c (revision 26674) @@ -126,5 +126,5 @@ int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits @@ -219,5 +219,5 @@ bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits Index: branches/eam_branches/20091201/ppStack/src/ppStackSetup.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackSetup.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackSetup.c (revision 26674) @@ -22,5 +22,7 @@ psAssert(recipe, "We've thrown an error on this before."); - options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis? + options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points? + + options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry? options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images? Index: branches/eam_branches/20091201/ppStack/src/ppStackSources.c =================================================================== --- branches/eam_branches/20091201/ppStack/src/ppStackSources.c (revision 26671) +++ branches/eam_branches/20091201/ppStack/src/ppStackSources.c (revision 26674) @@ -61,13 +61,13 @@ PS_ASSERT_PTR_NON_NULL(config, false); - if (!options->matchZPs) { - int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs - options->norm = psVectorAlloc(num, PS_TYPE_F32); - psVectorInit (options->norm, 0.0); - - // XXX do I need to set this? - // options->sumExposure = sumExpTime; - - return true; + if (!options->matchZPs && !options->photometry) { + int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + options->norm = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit (options->norm, 0.0); + + // XXX do I need to set this? + // options->sumExposure = sumExpTime; + + return true; } @@ -191,11 +191,68 @@ #endif - psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, - iter2, starRej2, starSys2, starLimit, - transIter, transRej, transThresh); // Transparencies for each image - if (!trans) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); - return false; - } + if (options->matchZPs) { + psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + iter2, starRej2, starSys2, starLimit, + transIter, transRej, transThresh); // Transparencies per image + if (!trans) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); + return false; + } + for (int i = 0; i < trans->n; i++) { + if (!isfinite(trans->data.F32[i])) { + inputMask->data.U8[i] |= PPSTACK_MASK_CAL; + } + } + // M = m + c0 + c1 * airmass - 2.5log(t) + transparency + // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 + // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 + // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) + // We don't need to know the magnitude zero point for the filter, since it cancels out + if (options->matchZPs) { + options->norm = psVectorAlloc(num, PS_TYPE_F32); + for (int i = 0; i < num; i++) { + if (!isfinite(trans->data.F32[i])) { + continue; + } + psArray *sources = sourceLists->data[i]; // Sources of interest + float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i]; + options->norm->data.F32[i] = magCorr; + psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", + i, magCorr); + + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source) { + continue; + } + source->psfMag += magCorr; + } + } + } + +#ifdef TESTING + dumpMatches("source_mags.dat", num, matches, zp, trans); +#endif + psFree(trans); + +#ifdef TESTING + // Double check: all transparencies should be zero + { + psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches + if (!matches) { + psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); + psFree(zp); + return false; + } + psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, + transThresh, starRej, starSys); + for (int i = 0; i < num; i++) { + fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); + } + psFree(trans); + } +#endif + } + #if 0 @@ -216,84 +273,28 @@ #endif -#ifdef TESTING - dumpMatches("source_mags.dat", num, matches, zp, trans); -#endif - - for (int i = 0; i < trans->n; i++) { - if (!isfinite(trans->data.F32[i])) { - inputMask->data.U8[i] |= PPSTACK_MASK_CAL; - } - } - - // Save best matches SOMEWHERE for future photometry - // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack - pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); - psArray *sourcesBest = psArrayAllocEmpty(matches->n); - - // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible - int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required - for (int i = 0; i < matches->n; i++) { - pmSourceMatch *match = matches->data[i]; // Match of interest - if (match->num < minMatches) { - continue; - } - - // We need to grab a single instance of this source: just take the first available - int image = match->image->data.S32[0]; // Index of image - int index = match->index->data.S32[0]; // Index of source within image - psArray *sources = sourceLists->data[image]; // Sources for image - pmSource *source = sources->data[index]; // Source of interest - - psArrayAdd(sourcesBest, sourcesBest->n, source); - } - psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, - "psphot sources", sourcesBest); - psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); - psFree(sourcesBest); + if (options->photometry) { + // Save best matches for future photometry + psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n); + + // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible + int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + if (match->num < minMatches) { + continue; + } + + // We need to grab a single instance of this source: just take the first available + int image = match->image->data.S32[0]; // Index of image + int index = match->index->data.S32[0]; // Index of source within image + psArray *sources = sourceLists->data[image]; // Sources for image + pmSource *source = sources->data[index]; // Source of interest + + psArrayAdd(sourcesBest, sourcesBest->n, source); + } + psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); + } psFree(matches); - - // M = m + c0 + c1 * airmass - 2.5log(t) + transparency - // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 - // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 - // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) - // We don't need to know the magnitude zero point for the filter, since it cancels out - options->norm = psVectorAlloc(num, PS_TYPE_F32); - for (int i = 0; i < num; i++) { - if (!isfinite(trans->data.F32[i])) { - continue; - } - psArray *sources = sourceLists->data[i]; // Sources of interest - float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i]; - options->norm->data.F32[i] = magCorr; - psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr); - - for (int j = 0; j < sources->n; j++) { - pmSource *source = sources->data[j]; // Source of interest - if (!source) { - continue; - } - source->psfMag += magCorr; - } - } - psFree(trans); - -#ifdef TESTING - // Double check: all transparencies should be zero - { - psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches - if (!matches) { - psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); - psFree(zp); - return false; - } - psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, - transThresh, starRej, starSys); - for (int i = 0; i < num; i++) { - fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); - } - psFree(trans); - } -#endif return true;