Index: branches/pap/psphot/src/psphotApResid.c =================================================================== --- branches/pap/psphot/src/psphotApResid.c (revision 28179) +++ branches/pap/psphot/src/psphotApResid.c (revision 28484) @@ -22,9 +22,9 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image - if (!psphotApResidReadout (config, view, filerule, i, recipe)) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotApResidReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -56,6 +56,6 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid"); + return true; } @@ -66,5 +66,5 @@ int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads if (!status) { - nThreads = 0; + nThreads = 0; } @@ -128,48 +128,46 @@ for (int i = 0; i < cellGroups->n; i++) { - psArray *cells = cellGroups->data[i]; - - for (int j = 0; j < cells->n; j++) { - - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS"); - - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nskip - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return false; - } - psFree(job); - } - - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[5]; - Nskip += scalar->data.S32; - scalar = job->args->data[6]; - Nfail += scalar->data.S32; - } - psFree(job); - } + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nskip + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Nskip += scalar->data.S32; + scalar = job->args->data[6]; + Nfail += scalar->data.S32; + } + psFree(job); + } } @@ -184,5 +182,5 @@ Npsf = 0; -# ifdef DEBUG +# ifdef DEBUG FILE *f = fopen ("apresid.dat", "w"); psAssert (f, "failed open"); @@ -199,14 +197,14 @@ if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED"); - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY"); - if (source->mode & PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT"); - + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED"); + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY"); + if (source->mode & PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT"); + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { continue; } - // XXX make this user-configurable? - if (source->errMag > 0.01) continue; + // XXX make this user-configurable? + if (source->errMag > 0.01) continue; // aperture residual for this source @@ -221,17 +219,17 @@ # ifdef DEBUG - fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", - source->peak->xf, source->peak->yf, - source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], - source->psfMag, source->apMag, source->errMag, - source->modelPSF->params->data.F32[PM_PAR_I0], - source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], - source->modelPSF->params->data.F32[PM_PAR_7]); + fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", + source->peak->xf, source->peak->yf, + source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], + source->psfMag, source->apMag, source->errMag, + source->modelPSF->params->data.F32[PM_PAR_I0], + source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], + source->modelPSF->params->data.F32[PM_PAR_7]); # endif - if (!isfinite(source->psfMag)) psAbort ("nan in psfMag"); - if (!isfinite(source->errMag)) psAbort ("nan in errMag"); - if (!isfinite(source->apMag)) psAbort ("nan in apMag"); - if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos"); - if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos"); + if (!isfinite(source->psfMag)) psAbort ("nan in psfMag"); + if (!isfinite(source->errMag)) psAbort ("nan in errMag"); + if (!isfinite(source->apMag)) psAbort ("nan in apMag"); + if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos"); + if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos"); psVectorAppend (mag, source->psfMag); @@ -253,9 +251,9 @@ if (Npsf < APTREND_NSTAR_MIN) { psWarning("Only %d valid aperture residual sources (need %d), giving up", Npsf, APTREND_NSTAR_MIN); - goto escape; - } - - // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different - // definition for 'order' (order_MAP = order_POLY + 1). in addition, we have a + goto escape; + } + + // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different + // definition for 'order' (order_MAP = order_POLY + 1). in addition, we have a // user-specified MAX order, which we should respect, regardless of the mode @@ -270,6 +268,6 @@ pmTrend2DMode mode = PM_TREND_MAP; if (mode == PM_TREND_MAP) { - MaxOrderForStars ++; - } + MaxOrderForStars ++; + } APTREND_ORDER_MAX = PS_MIN (APTREND_ORDER_MAX, MaxOrderForStars); @@ -283,37 +281,37 @@ int NY = readout->image->numRows; for (int i = 1; i <= APTREND_ORDER_MAX; i++) { - - int Nx, Ny; - if (NX > NY) { - Nx = i; - Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5)); - } else { - Ny = i; - Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5)); - } - - float errorFloor; + + int Nx, Ny; + if (NX > NY) { + Nx = i; + Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5)); + } else { + Ny = i; + Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5)); + } + + float errorFloor; pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag); - if (!apTrend) { - continue; - } - - // apply ApTrend results - // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; - // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; - // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center - // if (!isfinite(ApResid)) psAbort("nan apresid @ center"); + if (!apTrend) { + continue; + } + + // apply ApTrend results + // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5; + // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5; + // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center + // if (!isfinite(ApResid)) psAbort("nan apresid @ center"); // store the minimum errorFloor and best ApTrend to keep if (errorFloor < errorFloorMin) { errorFloorMin = errorFloor; - psFree (psf->ApTrend); - psf->ApTrend = psMemIncrRefCounter(apTrend); - } - psFree (apTrend); + psFree (psf->ApTrend); + psf->ApTrend = psMemIncrRefCounter(apTrend); + } + psFree (apTrend); } if (psf->ApTrend == NULL) { psWarning("Failed to find a valid aperture residual value"); - goto escape; + goto escape; } @@ -383,11 +381,11 @@ if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) { psWarning("Failed to fit trend for %d x %d map", Nx, Ny); - psFree (apTrend); - return NULL; + psFree (apTrend); + return NULL; } if (apTrend->mode == PM_TREND_MAP) { - // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST: - psImageMapRepair (apTrend->map->map); - // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST: + // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST: + psImageMapRepair (apTrend->map->map); + // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST: } @@ -400,6 +398,6 @@ if (!isfinite(*apResidSysErr)) { psWarning("Failed to find systematic error for %d x %d map", Nx, Ny); - psFree (apTrend); - return NULL; + psFree (apTrend); + return NULL; } @@ -408,6 +406,6 @@ if (psTraceGetLevel("psphot") >= 4) { - char filename[64]; - snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny); + char filename[64]; + snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny); FILE *dumpFile = fopen (filename, "w"); for (int i = 0; i < xPos->n; i++) { @@ -457,28 +455,28 @@ } - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - - bool status = pmSourceMagnitudes (source, psf, photMode, maskVal); - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); + + bool status = pmSourceMagnitudes (source, psf, photMode, maskVal); + + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - if (!status) { - Nskip ++; - psTrace ("psphot", 3, "skip : bad source mag"); - continue; - } - - if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { - Nfail ++; - psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); - continue; - } - source->mode |= PM_SOURCE_MODE_AP_MAGS; + if (!status) { + Nskip ++; + psTrace ("psphot", 3, "skip : bad source mag"); + continue; + } + + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { + Nfail ++; + psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); + continue; + } + source->mode |= PM_SOURCE_MODE_AP_MAGS; } Index: branches/pap/psphot/src/psphotBlendFit.c =================================================================== --- branches/pap/psphot/src/psphotBlendFit.c (revision 28179) +++ branches/pap/psphot/src/psphotBlendFit.c (revision 28484) @@ -15,8 +15,8 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) { + if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -48,6 +48,6 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit"); + return true; } @@ -87,6 +87,6 @@ sources = psArraySort (sources, pmSourceSortBySN); if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; } @@ -103,26 +103,24 @@ for (int j = 0; j < cells->n; j++) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); - psArray *newSources = psArrayAllocEmpty(16); - - psArrayAdd(job->args, 1, readout); - psArrayAdd(job->args, 1, recipe); - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - psArrayAdd(job->args, 1, newSources); // return for new sources - psFree (newSources); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return NULL; - } - psFree(job); + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); + psArray *newSources = psArrayAllocEmpty(16); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, recipe); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, newSources); // return for new sources + psFree (newSources); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } # if (0) @@ -152,33 +150,33 @@ } - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[5]; - Nfit += scalar->data.S32; - scalar = job->args->data[6]; - Npsf += scalar->data.S32; - scalar = job->args->data[7]; - Next += scalar->data.S32; - scalar = job->args->data[8]; - Nfail += scalar->data.S32; - - // add these back onto sources - psArray *newSources = job->args->data[4]; - for (int j = 0; j < newSources->n; j++) { - psArrayAdd (sources, 16, newSources->data[j]); - } - } - psFree(job); + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Nfit += scalar->data.S32; + scalar = job->args->data[6]; + Npsf += scalar->data.S32; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nfail += scalar->data.S32; + + // add these back onto sources + psArray *newSources = job->args->data[4]; + for (int j = 0; j < newSources->n; j++) { + psArrayAdd (sources, 16, newSources->data[j]); + } + } + psFree(job); } } @@ -278,5 +276,5 @@ psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); Next ++; - source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; continue; } @@ -286,5 +284,5 @@ psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); Npsf ++; - source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; continue; } Index: branches/pap/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/pap/psphot/src/psphotFitSourcesLinear.c (revision 28179) +++ branches/pap/psphot/src/psphotFitSourcesLinear.c (revision 28484) @@ -27,24 +27,24 @@ for (int i = 0; i < num; i++) { - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - psAssert (psf, "missing psf?"); - - if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) { psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -59,6 +59,6 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); + return true; } @@ -150,4 +150,5 @@ if (fitSources->n == 0) { + psFree(fitSources); return true; } @@ -170,10 +171,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor); + float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -183,5 +184,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseVectorElement (sparse, i, f); @@ -189,11 +190,11 @@ switch (SKY_FIT_ORDER) { case 1: - f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 1, f); - f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 2, f); case 0: - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseBorderElementB (border, i, 0, f); break; @@ -215,5 +216,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal); psSparseMatrixElement (sparse, j, i, f); } Index: branches/pap/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/pap/psphot/src/psphotFitSourcesLinearStack.c (revision 28179) +++ branches/pap/psphot/src/psphotFitSourcesLinearStack.c (revision 28484) @@ -43,25 +43,25 @@ for (int i = 0; i < objects->n; i++) { pmPhotObj *object = objects->data[i]; - if (!object) continue; - if (!object->sources) continue; + if (!object) continue; + if (!object->sources) continue; - // XXX check an element of the group to see if we should use it - // if (!object->flags & PM_PHOT_OBJ_BAD) continue; + // XXX check an element of the group to see if we should use it + // if (!object->flags & PM_PHOT_OBJ_BAD) continue; - for (int j = 0; j < object->sources->n; j++) { - pmSource *source = object->sources->data[j]; - if (!source) continue; + for (int j = 0; j < object->sources->n; j++) { + pmSource *source = object->sources->data[j]; + if (!source) continue; - // turn this bit off and turn it on again if we keep this source - source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + // turn this bit off and turn it on again if we keep this source + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; - // generate model for sources without, or skip if we can't - if (!source->modelFlux) { + // generate model for sources without, or skip if we can't + if (!source->modelFlux) { if (!pmSourceCacheModel (source, maskVal)) continue; - } + } - source->mode |= PM_SOURCE_MODE_LINEAR_FIT; - psArrayAdd (fitSources, 100, source); - } + source->mode |= PM_SOURCE_MODE_LINEAR_FIT; + psArrayAdd (fitSources, 100, source); + } } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); @@ -85,10 +85,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR); + float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -97,5 +97,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseVectorElement (sparse, i, f); @@ -104,6 +104,6 @@ pmSource *SRCj = fitSources->data[j]; - // we only need to generate dot terms for source on the same image - if (SRCj->imageID != SRCi->imageID) { continue; } + // we only need to generate dot terms for source on the same image + if (SRCj->imageID != SRCi->imageID) { continue; } // skip over disjoint source images, break after last possible overlap @@ -114,5 +114,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal); psSparseMatrixElement (sparse, j, i, f); } Index: branches/pap/psphot/src/psphotGuessModels.c =================================================================== --- branches/pap/psphot/src/psphotGuessModels.c (revision 28179) +++ branches/pap/psphot/src/psphotGuessModels.c (revision 28484) @@ -21,5 +21,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image + if (i == chisqNum) continue; // skip chisq image if (!psphotGuessModelsReadout (config, view, filerule, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for %s entry %d", filerule, i); @@ -106,8 +106,6 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); return false; } - psFree(job); } Index: branches/pap/psphot/src/psphotImageLoop.c =================================================================== --- branches/pap/psphot/src/psphotImageLoop.c (revision 28179) +++ branches/pap/psphot/src/psphotImageLoop.c (revision 28484) @@ -46,4 +46,21 @@ if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); + // Read WCS if easy. + // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. + { + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip + pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + // try to load other supporting data (PSF, SRC, etc). // do not re-load the following three files @@ -67,19 +84,21 @@ // Update the header - pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); - if (hdu && hdu != lastHDU) { - psphotVersionHeaderFull(hdu->header); - lastHDU = hdu; + { + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } } - // if an external mask is supplied, ensure that NAN pixels are also masked - if (readout->mask) { - psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels - if (!pmReadoutMaskNonfinite(readout, maskSat)) { - psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); - psFree(view); - return false; - } - } + // if an external mask is supplied, ensure that NAN pixels are also masked + if (readout->mask) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels + if (!pmReadoutMaskNonfinite(readout, maskSat)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); + psFree(view); + return false; + } + } // run the actual photometry analysis on this chip/cell/readout @@ -91,14 +110,14 @@ } - // drop all versions of the internal files - status = true; - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); - if (!status) { - psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); - psFree (view); - return false; - } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } } // save output which is saved at the chip level Index: branches/pap/psphot/src/psphotMagnitudes.c =================================================================== --- branches/pap/psphot/src/psphotMagnitudes.c (revision 28179) +++ branches/pap/psphot/src/psphotMagnitudes.c (revision 28484) @@ -18,26 +18,26 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - psAssert (psf, "missing psf?"); - - if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) { + if (i == chisqNum) continue; // skip chisq image + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -50,8 +50,8 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes"); - return true; - } - + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes"); + return true; + } + psTimerStart ("psphot.mags"); @@ -122,8 +122,6 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); return false; } - psFree(job); # if (0) @@ -184,13 +182,13 @@ } - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); + // clear the mask bit and set the circular mask pixels + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal if (status && isfinite(source->apMag)) Nap ++; - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + // clear the mask bit + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // re-subtract the object, leave local sky @@ -273,9 +271,6 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); return false; } - psFree(job); - } Index: branches/pap/psphot/src/psphotRadiusChecks.c =================================================================== --- branches/pap/psphot/src/psphotRadiusChecks.c (revision 28179) +++ branches/pap/psphot/src/psphotRadiusChecks.c (revision 28484) @@ -4,7 +4,7 @@ static float PSF_FIT_NSIGMA; static float PSF_FIT_PADDING; -static float PSF_APERTURE = 0; // radius to use in PSF aperture mags -static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, - // and a per-object radius is calculated) +static float PSF_APERTURE = 0; // radius to use in PSF aperture mags +static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, + // and a per-object radius is calculated) bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { @@ -17,10 +17,30 @@ PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); if (!status) { - PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); + PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); } PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); if (!status) { - PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + } + + // The PSF_FIT_RADIUS and PSF_APERTURE may not be set if the PSF was loaded and not chosen + + if (PSF_FIT_RADIUS == 0.0) { + float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE"); + PSF_FIT_RADIUS = (int)(fitScale*gaussSigma); + } + + if (PSF_APERTURE == 0.0) { + float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); + PSF_APERTURE = (int)(apScale*gaussSigma); } @@ -38,18 +58,18 @@ // set the fit radius based on the object flux limit and the model float radiusFit = PSF_FIT_RADIUS; - if (radiusFit <= 0) { // use fixed radius - if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); - } else { - radiusFit = model->modelRadius(model->params, 1.0); - } - model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + if (radiusFit <= 0) { // use fixed radius + if (moments == NULL) { + radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + } else { + radiusFit = model->modelRadius(model->params, 1.0); + } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); } else { - model->fitRadius = radiusFit; + model->fitRadius = radiusFit; } if (isnan(model->fitRadius)) psAbort("error in radius"); - + if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->fitRadius *= 2; + model->fitRadius *= 2; } @@ -73,13 +93,13 @@ // set the fit radius based on the object flux limit and the model float radiusFit = PSF_FIT_RADIUS; - if (radiusFit <= 0) { // use fixed radius - if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); - } else { - radiusFit = model->modelRadius(model->params, 1.0); - } - model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); + if (radiusFit <= 0) { // use fixed radius + if (moments == NULL) { + radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + } else { + radiusFit = model->modelRadius(model->params, 1.0); + } + model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); } else { - model->fitRadius = radiusFit; + model->fitRadius = radiusFit; } if (isnan(model->fitRadius)) psAbort("error in radius"); @@ -89,5 +109,5 @@ if (source->mode & PM_SOURCE_MODE_SATSTAR) { - model->fitRadius *= 2; + model->fitRadius *= 2; } @@ -134,12 +154,12 @@ float radius = 0.0; for (int j = 0; j < footprint->spans->n; j++) { - pmSpan *span = footprint->spans->data[j]; - - float dY = span->y - peak->yf; - float dX0 = span->x0 - peak->xf; - float dX1 = span->x1 - peak->xf; - - radius = PS_MAX (radius, hypot(dY, dX0)); - radius = PS_MAX (radius, hypot(dY, dX1)); + pmSpan *span = footprint->spans->data[j]; + + float dY = span->y - peak->yf; + float dX0 = span->x0 - peak->xf; + float dX1 = span->x1 - peak->xf; + + radius = PS_MAX (radius, hypot(dY, dX0)); + radius = PS_MAX (radius, hypot(dY, dX1)); } Index: branches/pap/psphot/src/psphotReadout.c =================================================================== --- branches/pap/psphot/src/psphotReadout.c (revision 28179) +++ branches/pap/psphot/src/psphotReadout.c (revision 28484) @@ -45,4 +45,5 @@ return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } + if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); @@ -93,4 +94,5 @@ return psphotReadoutCleanup (config, view, "PSPHOT.INPUT"); } + // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) { // pass 1 @@ -208,6 +210,8 @@ // calculate source magnitudes - psphotMagnitudes(config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources) - + if (!psphotMagnitudes(config, view, "PSPHOT.INPUT")) { // pass 1 (detections->allSources) + psErrorStackPrint(stderr, "Unable to do magnitudes."); + psErrorClear(); + } if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) { // pass 1 psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); @@ -219,9 +223,19 @@ // replace background in residual image - psphotSkyReplace (config, view, "PSPHOT.INPUT"); // pass 1 - + if (!psphotSkyReplace (config, view, "PSPHOT.INPUT")) { // pass 1 + psErrorStackPrint(stderr, "Unable to replace sky"); + psErrorClear(); + +/* psLogMsg("psphot", 3, "failed on psphotSkyReplace"); */ +/* return(psphotReadoutCleanup(config, view, "PSPHOT.INPUT")); */ + } // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // pass 1 - + if (!psphotSourceFreePixels (config, view, "PSPHOT.INPUT")) { // pass 1 + psErrorStackPrint(stderr, "Unable to free source pixels"); + psErrorClear(); + +/* psLogMsg ("psphot", 3, "failed on psphotSourceFreePixels"); */ +/* return(psphotReadoutCleanup(config, view, "PSPHOT.INPUT")); */ + } // create the exported-metadata and free local data return psphotReadoutCleanup(config, view, "PSPHOT.INPUT"); Index: branches/pap/psphot/src/psphotRoughClass.c =================================================================== --- branches/pap/psphot/src/psphotRoughClass.c (revision 28179) +++ branches/pap/psphot/src/psphotRoughClass.c (revision 28484) @@ -25,9 +25,9 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (i == chisqNum) continue; // skip chisq image - if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i); - return false; - } + return false; + } } return true; @@ -50,5 +50,5 @@ bool havePSF = false; if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { - havePSF = true; + havePSF = true; } @@ -60,6 +60,6 @@ if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); - return true; + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); + return true; } @@ -78,4 +78,17 @@ psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + + // If in doubt, it's a PSF + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source || !source->peak) { + continue; + } + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + source->type = PM_SOURCE_TYPE_STAR; + } psFree (region); continue; @@ -124,22 +137,22 @@ // XXX why not save the psfClump as a PTR? - float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); - float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); + float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); - float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); - if (!status) { - PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); - psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); - } - float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); - if (!status) { - MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); - psAssert (status, "missing MOMENTS_SX_MAX"); - } - float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); - if (!status) { - MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); - psAssert (status, "missing MOMENTS_SY_MAX"); - } + float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); + if (!status) { + PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); + psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); + } + float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); + if (!status) { + MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); + psAssert (status, "missing MOMENTS_SX_MAX"); + } + float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); + if (!status) { + MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); + psAssert (status, "missing MOMENTS_SY_MAX"); + } psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); Index: branches/pap/psphot/src/psphotSourceStats.c =================================================================== --- branches/pap/psphot/src/psphotSourceStats.c (revision 28179) +++ branches/pap/psphot/src/psphotSourceStats.c (revision 28484) @@ -3,5 +3,5 @@ // convert detections to sources and measure their basic properties (moments, local sky, sky // variance) Note: this function only generates sources for the new peaks (peak->assigned). -// The new sources are added to any existing sources on detections->newSources. The sources +// The new sources are added to any existing sources on detections->newSources. The sources // on detections->allSources are ignored. bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow) @@ -22,5 +22,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - // if (i == chisqNum) continue; // skip chisq image + // if (i == chisqNum) continue; // skip chisq image if (!psphotSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) { psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i); @@ -91,5 +91,5 @@ // generate the array of sources, define the associated pixel if (!detections->newSources) { - detections->newSources = psArrayAllocEmpty (peaks->n); + detections->newSources = psArrayAllocEmpty (peaks->n); } sources = detections->newSources; @@ -107,5 +107,5 @@ // create a new source pmSource *source = pmSourceAlloc(); - source->imageID = index; + source->imageID = index; // add the peak @@ -184,9 +184,7 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); - psFree (job); psFree(detections->newSources); return false; } - psFree(job); } @@ -194,5 +192,5 @@ if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); - psFree(detections->newSources); + psFree(detections->newSources); return false; } @@ -315,8 +313,6 @@ if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); - psFree (job); return NULL; } - psFree(job); } @@ -380,6 +376,6 @@ pmSource *source = sources->data[i]; - if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue; - source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue; + source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; // skip faint sources for moments measurement @@ -493,4 +489,6 @@ psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); +#if 0 + // Modifying clump parameters without restoring! psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", 1); psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000"); @@ -504,6 +502,6 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); - psphotVisualPlotMoments (recipe, analysis, sources); +#endif Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; Index: branches/pap/psphot/src/psphotTest.c =================================================================== --- branches/pap/psphot/src/psphotTest.c (revision 28179) +++ branches/pap/psphot/src/psphotTest.c (revision 28484) @@ -25,5 +25,5 @@ bool FillImage_Threaded (psThreadJob *job) { - + psImage *image = job->args->data[0]; int xs = PS_SCALAR_VALUE(job->args->data[1],S32); @@ -36,7 +36,7 @@ psRegion region = psRegionSet (xs, xs + dx, ys, ys + dy); for (int i = 0; i < 100; i++) { - psImage *subset = psImageSubset (image, region); - psImageInit (subset, value + i); - psFree (subset); + psImage *subset = psImageSubset (image, region); + psImageInit (subset, value + i); + psFree (subset); } return true; @@ -46,6 +46,6 @@ if (argc != 3) { - fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n"); - exit (2); + fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n"); + exit (2); } @@ -62,25 +62,23 @@ for (int ix = 0; ix < 1000; ix += 100) { - for (int iy = 0; iy < 1000; iy += 100) { + for (int iy = 0; iy < 1000; iy += 100) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE"); + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE"); - psArrayAdd(job->args, 1, image); - PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32); + psArrayAdd(job->args, 1, image); + PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32); - // FillImage (image, ix, iy, 100, 100, ix + iy); + // FillImage (image, ix, iy, 100, 100, ix + iy); - if (!psThreadJobAddPending(job)) { - fprintf (stderr, "failure to run FillImage(1)"); - psFree (job); - exit (1); - } - psFree(job); - } + if (!psThreadJobAddPending(job)) { + fprintf (stderr, "failure to run FillImage(1)"); + exit (1); + } + } } @@ -88,6 +86,6 @@ // wait for the threads to finish and manage results if (!psThreadPoolWait (true)) { - fprintf (stderr, "failure to run FillImage (2)"); - exit (1); + fprintf (stderr, "failure to run FillImage (2)"); + exit (1); }