Index: /trunk/ippconfig/recipes/ppStack.config =================================================================== --- /trunk/ippconfig/recipes/ppStack.config (revision 23211) +++ /trunk/ippconfig/recipes/ppStack.config (revision 23212) @@ -32,12 +32,16 @@ ZP.RADIUS F32 1.0 # Radius (pixels) for matching sources -ZP.ITER S32 1000 # Maximum iterations for zero point +ZP.ITER.1 S32 5 # Iterations for zero point; pass 1 +ZP.ITER.2 S32 1000 # Iterations for zero point; pass 2 ZP.TOL F32 1.0e-6 # Tolerance for zero point iterations ZP.TRANS.ITER S32 2 # Iterations for transparency determination ZP.TRANS.REJ F32 3.0 # Rejection threshold for transparency determination ZP.TRANS.THRESH F32 1.0 # Threshold for transparency determination -ZP.STAR.REJ F32 3.0 # Rejection threshold for stars +ZP.STAR.REJ.1 F32 20.0 # Rejection threshold for stars; pass 1 +ZP.STAR.REJ.2 F32 3.0 # Rejection threshold for stars; pass 2 ZP.STAR.LIMIT F32 1.0e-2 # Limit on star rejection fraction for successful iteration -ZP.STAR.SYS F32 0.05 # Estimated systematic error +ZP.STAR.SYS.1 F32 0.10 # Estimated systematic error; pass 1 +ZP.STAR.SYS.2 F32 0.05 # Estimated systematic error; pass 2 +ZP.MATCH F32 0.3 # Fraction of images to match for good star ZP.AIRMASS METADATA # Airmass terms by filter g F32 0.0 Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 23211) +++ /trunk/ppStack/src/ppStackArguments.c (revision 23212) @@ -168,12 +168,15 @@ "Play safe with small numbers of pixels to combine?", false); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-radius", 0, "Radius (pixels) for matching sources", NAN); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter", 0, "Maximum iterations for zero point", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter-1", 0, "Maximum iterations for zero point; pass 1", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter-2", 0, "Maximum iterations for zero point; pass 2", 0); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-tol", 0, "Tolerance for zero point iterations", NAN); psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-trans-iter", 0, "Iterations for transparency determination", 0); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-rej", 0, "Rejection threshold for transparency determination", NAN); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-thresh", 0, "Threshold for transparency determination", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej", 0, "Rejection threshold for stars", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej-1", 0, "Rejection threshold for stars; pass 1", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej-2", 0, "Rejection threshold for stars; pass 2", NAN); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-limit", 0, "Limit on star rejection fraction for successful iteration", NAN); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys", 0, "Estimated systematic error", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-1", 0, "Estimated systematic error; pass 1", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-2", 0, "Estimated systematic error; pass 2", NAN); psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL); @@ -242,12 +245,15 @@ VALUE_ARG_RECIPE_FLOAT("-zp-radius", "ZP.RADIUS", F32); - VALUE_ARG_RECIPE_INT("-zp-iter", "ZP.ITER", S32, 0); + VALUE_ARG_RECIPE_INT("-zp-iter-1", "ZP.ITER.1", S32, 0); + VALUE_ARG_RECIPE_INT("-zp-iter-2", "ZP.ITER.2", S32, 0); VALUE_ARG_RECIPE_FLOAT("-zp-tol", "ZP.TOL", F32); VALUE_ARG_RECIPE_INT("-zp-trans-iter", "ZP.TRANS.ITER", S32, 0); VALUE_ARG_RECIPE_FLOAT("-zp-trans-rej", "ZP.TRANS.REJ", F32); VALUE_ARG_RECIPE_FLOAT("-zp-trans-thresh", "ZP.TRANS.THRESH", F32); - VALUE_ARG_RECIPE_FLOAT("-zp-star-rej", "ZP.STAR.REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-zp-star-rej-1", "ZP.STAR.REJ.1", F32); + VALUE_ARG_RECIPE_FLOAT("-zp-star-rej-2", "ZP.STAR.REJ.2", F32); VALUE_ARG_RECIPE_FLOAT("-zp-star-limit", "ZP.STAR.LIMIT", F32); - VALUE_ARG_RECIPE_FLOAT("-zp-star-sys", "ZP.STAR.SYS", F32); + VALUE_ARG_RECIPE_FLOAT("-zp-star-sys-1", "ZP.STAR.SYS.1", F32); + VALUE_ARG_RECIPE_FLOAT("-zp-star-sys-2", "ZP.STAR.SYS.2", F32); VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0); Index: /trunk/ppStack/src/ppStackSources.c =================================================================== --- /trunk/ppStack/src/ppStackSources.c (revision 23211) +++ /trunk/ppStack/src/ppStackSources.c (revision 23212) @@ -81,12 +81,19 @@ float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources - int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations + int iter1 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.1"); // Maximum iterations for pass 1 + int iter2 = psMetadataLookupS32(NULL, recipe, "ZP.ITER.2"); // Maximum iterations for pass 2 float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency - float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars + + float starRej1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.1"); // Rejection threshold for stars + float starSys1 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.1"); // Estimated systematic error + float starRej2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ.2"); // Rejection threshold for stars + float starSys2 = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS.2"); // Estimated systematic error + float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction - float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error + + float fracMatch = psMetadataLookupF32(NULL, recipe, "ZP.MATCH"); // Fraction of images to match for star psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms @@ -164,6 +171,7 @@ #endif - psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, - transThresh, starRej, starSys); // Transparencies for each image + 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"); @@ -184,8 +192,8 @@ // 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 (100); - - // XXX something of a hack: require at least 2 detections or 1/2 of the max possible - int minMatches = PS_MAX(2, 0.3*num);// Minimum number of matches required + 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 @@ -200,5 +208,5 @@ pmSource *source = sources->data[index]; // Source of interest - psArrayAdd(sourcesBest, 100, source); + psArrayAdd(sourcesBest, sourcesBest->n, source); } psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, Index: /trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.c (revision 23211) +++ /trunk/psModules/src/objects/pmSourceMatch.c (revision 23212) @@ -524,12 +524,15 @@ psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches const psVector *zp, // Zero points for each image (including airmass term) - int maxIter, // Maximum number of iterations float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + float sys1, // Systematic error in measurements for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float sys2, // Systematic error in measurements for pass 2 float rejLimit, // Limit on rejection between iterations int transIter, // Clipping iterations for transparency float transClip, // Clipping level for transparency - float photoLevel, // Level at which we declare image is photometric - float starClip, // Clipping for stars - float sysErr // Systematic error in measurements + float photoLevel // Level at which we declare image is photometric ) { @@ -539,5 +542,6 @@ PS_ASSERT_FLOAT_LARGER_THAN(transClip, 0.0, NULL); - sysErr *= sysErr; + sys1 *= sys1; + sys2 *= sys2; int numImages = zp->n; // Number of images @@ -552,11 +556,12 @@ float chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, - photo, sysErr); // chi^2 for solution + photo, sys1); // chi^2 for solution psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2); float lastChi2 = INFINITY; // chi^2 on last iteration float fracRej = INFINITY; // Fraction of measurements rejected - // in the first passes, the transparencies are not well deteremined: use high systematic error and the rejection thresholds - for (int i = 0; i < 5; i++) { + // In the first passes, the transparencies are not well deteremined: use high systematic error and + // rejection thresholds + for (int i = 0; i < iter1; i++) { // Identify photometric nights @@ -572,14 +577,12 @@ psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages); - // XXX use 20 sigma rejection and 0.1 mag systematic error (move these to the recipe) - fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, 20.0, PS_SQR(0.1)); + fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1); psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100); - // XXX use 0.05 mag systematic error (move these to the recipe) - chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, PS_SQR(0.1)); + chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1); psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); } - for (int i = 0; i < maxIter && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) { + for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) { lastChi2 = chi2; @@ -594,11 +597,11 @@ return NULL; } - psTrace("psModules.objects", 3, "Determined %d/%d are photometric", numPhoto, numImages); - - fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, starClip, sysErr); - psTrace("psModules.objects", 3, "%f%% of measurements rejected", fracRej * 100); - - chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sysErr); - psTrace("psModules.objects", 1, "iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages); + + fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2); + psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100); + + chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2); + psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); } Index: /trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- /trunk/psModules/src/objects/pmSourceMatch.h (revision 23211) +++ /trunk/psModules/src/objects/pmSourceMatch.h (revision 23212) @@ -55,12 +55,15 @@ psVector *pmSourceMatchRelphot(const psArray *matches, // Array of matches const psVector *zp, // Zero points for each image (including airmass term) - int maxIter, // Maximum number of iterations float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + float sys1, // Systematic error in measurements for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float sys2, // Systematic error in measurements for pass 2 float rejLimit, // Limit on rejection between iterations int transIter, // Clipping iterations for transparency float transClip, // Clipping level for transparency - float photoLevel, // Level at which we declare image is photometric - float starClip, // Clipping for stars - float sysErr // Systematic error in measurements + float photoLevel // Level at which we declare image is photometric );