Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 6570) +++ /trunk/psphot/src/Makefile.am (revision 6571) @@ -12,7 +12,12 @@ psphot_SOURCES = \ + psphot.c \ + psphotModelGroupInit.c \ + psphotArguments.c \ + psphotParseCamera.c \ + testPSphotLoop.c \ + psphotMaskCell.c \ psphotReadout.c \ - psphotImageStats.c \ - psphotImageBackground.c \ + psphotImageMedian.c \ psphotFindPeaks.c \ psphotSourceStats.c \ @@ -24,31 +29,20 @@ psphotReplaceUnfit.c \ psphotApResid.c \ - psphotOutput.c \ - psphotLoadPixels.c \ - psphotModelGroupInit.c \ - psphotGrowthCurve.c \ - psphotCleanup.c \ + psphotMagnitudes.c \ + psphotSkyReplace.c \ psphotEvalPSF.c \ psphotEvalFLT.c \ - psphotSourceFits.c \ + psphotSourceFits.c + +junk = psphotGrowthCurve.c \ + psphotCleanup.c \ psphotSortBySN.c \ psphotDefinePixels.c \ - psphotMagnitudes.c \ psphotRadiusChecks.c \ - psphotImageMedian.c \ - psphotSkyReplace.c \ - psLine.c \ - psSparse.c \ - psPolynomialUtils.c \ psModulesUtils.c \ pmSourceContour.c \ pmModelFitSet.c \ - pmCellSetMask.c \ pmSourceFitSet.c \ - psphot.c \ - psphotArguments.c \ - psphotParseCamera.c \ - psphotModelTest.c \ - psphotImageLoop.c + psphotModelTest.c noinst_HEADERS = \ @@ -68,2 +62,6 @@ # psphotApplyPSF.c # psphotFitGalaxies.c +# psphotImageStats.c +# psphotImageBackground.c +# psphotOutput.c +# psphotLoadPixels.c Index: unk/psphot/src/pmCellSetMask.c =================================================================== --- /trunk/psphot/src/pmCellSetMask.c (revision 6570) +++ (revision ) @@ -1,31 +1,0 @@ -# include "psphot.h" - -// 2006.02.04 : no leaks -bool pmCellSetMask (pmCell *cell, psMetadata *recipe) { - - bool status; - - // mask the excluded outer pixels - // these coordinates refer to the parent image - // these bounds will saturate on the subimage - // negative upper bounds will subtract from the *subimage* - float XMIN = psMetadataLookupF32 (&status, recipe, "XMIN"); - float XMAX = psMetadataLookupF32 (&status, recipe, "XMAX"); - float YMIN = psMetadataLookupF32 (&status, recipe, "YMIN"); - float YMAX = psMetadataLookupF32 (&status, recipe, "YMAX"); - psRegion valid = psRegionSet (XMIN, XMAX, YMIN, YMAX); - - // Set the pixels - psArray *readouts = cell->readouts; // Array of readouts - for (int i = 0; i < readouts->n; i++) { - pmReadout *readout = readouts->data[i]; // The readout of interest - - // psRegionForImage creates a region for the subimage, not the parent - // XXX I propose that psRegionForImage treat the region as referring to the parent... - psRegion keep = psRegionForImage (readout->image, valid); - - // psImageKeepRegion assumes the region refers to the parent coordinates - psImageKeepRegion (readout->mask, keep, "OR", PSPHOT_MASK_INVALID); - } - return true; -} Index: unk/psphot/src/pmModelFitSet.c =================================================================== --- /trunk/psphot/src/pmModelFitSet.c (revision 6570) +++ (revision ) @@ -1,72 +1,0 @@ -# include "psphot.h" - -// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, -// nPar = nSrc*(nOnePar - 1) + 1 - -static pmModelFunc mFunc; -static int nPar; -static psVector *onePar; -static psVector *oneDeriv; - -bool pmModelFitSetInit (pmModelType type) { - - mFunc = pmModelFunc_GetFunction (type); - nPar = pmModelParameterCount (type); - - onePar = psVectorAlloc (nPar, PS_DATA_F32); - oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); - - return true; -} - -void pmModelFitSetClear (void) { - - psFree (onePar); - psFree (oneDeriv); - return; -} - -psF32 pmModelFitSet(psVector *deriv, - const psVector *params, - const psVector *x) -{ - - psF32 value; - psF32 model; - - psF32 *PAR = onePar->data.F32; - psF32 *dPAR = oneDeriv->data.F32; - - psF32 *pars = params->data.F32; - psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; - - int nSrc = (params->n - 1) / (nPar - 1); - - PAR[0] = model = pars[0]; - for (int i = 0; i < nSrc; i++) { - int nOff = i*nPar - i; - for (int n = 1; n < nPar; n++) { - PAR[n] = pars[nOff + n]; - } - if (deriv == NULL) { - value = mFunc (NULL, onePar, x); - } else { - value = mFunc (oneDeriv, onePar, x); - for (int n = 1; n < nPar; n++) { - dpars[nOff + n] = dPAR[n]; - } - } - model += value; - } - if (deriv != NULL) { - dpars[0] = dPAR[0]*2.0; - } - return (model); -} - -/* -i: 0 1 2 -n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 -i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 -*/ - Index: unk/psphot/src/pmSourceContour.c =================================================================== --- /trunk/psphot/src/pmSourceContour.c (revision 6570) +++ (revision ) @@ -1,213 +1,0 @@ -# include "psphot.h" - -/****************************************************************************** -findValue(source, level, row, col, dir): a private function which determines -the column coordinate of the model function which has the value "level". If -dir equals 0, then you loop leftwards from the peak pixel, otherwise, -rightwards. - -XXX: reverse order of row,col args? - -XXX: Input row/col are in image coords. - -XXX: The result is returned in image coords. -*****************************************************************************/ -# define LEFT false -# define RIGHT true - -// return the first coordinate at or below the threshold in the requested direction -static int findContourNeg( - psImage *image, - float threshold, - int x, - int y, - bool right) -{ - - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - - // We define variables incr and lastColumn so that we can use the same loop - // whether we are stepping leftwards, or rightwards. - - int incr; - int subCol; - int lastColumn; - if (right) { - incr = 1; - lastColumn = image->numCols - 1; - } else { - incr = -1; - lastColumn = 0; - } - - subCol = x; - - while (subCol != lastColumn) { - float value = image->data.F32[y][subCol]; - if (value <= threshold) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return (subCol); - } - subCol += incr; - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return (lastColumn); -} - -// return the last coordinate at or below the threshold in the requested direction -static int findContourPos( - psImage *image, - float threshold, - int x, - int y, - bool right, - int xEnd) -{ - - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - - // We define variables incr and lastColumn so that we can use the same loop - // whether we are stepping leftwards, or rightwards. - - int incr; - int subCol; - int lastColumn; - if (right) { - incr = 1; - lastColumn = PS_MIN (image->numCols - 1, xEnd); - } else { - incr = -1; - lastColumn = PS_MAX (0, xEnd); - } - - subCol = x; - while (subCol != lastColumn) { - float value = image->data.F32[y][subCol]; - if (value >= threshold) { - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - if (subCol == x) { - return (subCol); - } - return (subCol); - } - subCol += incr; - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return (lastColumn); -} - -/****************************************************************************** -new implementation of source contour function -*****************************************************************************/ -psArray *pmSourceContour_EAM (psImage *image, int xc, int yc, float threshold) { - - int xR, yR, x0, x1, x0s, x1s; - - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(image, false); - - int x = xc - image->col0; - int y = yc - image->row0; - - // the requested point must be within the contour - if (image->data.F32[y][x] < threshold) return NULL; - - // Ensure that the starting column is allowable. - if (x < 0) return NULL; - if (y < 0) return NULL; - if (x >= image->numCols) return NULL; - if (y >= image->numRows) return NULL; - - // Allocate data for x/y pairs. - psVector *xVec = psVectorAlloc(100, PS_TYPE_F32); - psVector *yVec = psVectorAlloc(100, PS_TYPE_F32); - - // First row: find the left and right end-points - int Npt = 0; - - x0 = findContourNeg (image, threshold, x, y, LEFT); - x1 = findContourNeg (image, threshold, x, y, RIGHT); - xVec->data.F32[Npt + 0] = image->col0 + x0; - xVec->data.F32[Npt + 1] = image->col0 + x1; - yVec->data.F32[Npt + 0] = image->row0 + y; - yVec->data.F32[Npt + 1] = image->row0 + y; - Npt += 2; - - x0s = x0; - x1s = x1; - - // look for contour outline above row - xR = x0s; - yR = y + 1; - while (yR < image->numRows) { - if (image->data.F32[yR][xR] < threshold) { - x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); - if (x0 == x1) { - // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); - goto pt1; - } - x1 = findContourNeg (image, threshold, x0, yR, RIGHT); - x0--; - } else { - x0 = findContourNeg (image, threshold, xR, yR, LEFT); - x1 = findContourNeg (image, threshold, xR, yR, RIGHT); - } - // fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); - - xVec->data.F32[Npt + 0] = image->col0 + x0; - xVec->data.F32[Npt + 1] = image->col0 + x1; - yVec->data.F32[Npt + 0] = image->row0 + yR; - yVec->data.F32[Npt + 1] = image->row0 + yR; - Npt += 2; - - if (Npt >= xVec->nalloc - 1) { - psVectorRealloc (xVec, xVec->nalloc + 100); - psVectorRealloc (yVec, yVec->nalloc + 100); - } - yR ++; - xR = x0; - } - -pt1: - // look for contour outline below row - xR = x0s; - x1 = x1s; - yR = y - 1; - while (yR >= 0) { - if (image->data.F32[yR][xR] < threshold) { - x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); - if (x0 == x1) { - // fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); - goto pt2; - } - x1 = findContourNeg (image, threshold, x0, yR, RIGHT); - x0--; - } else { - x0 = findContourNeg (image, threshold, xR, yR, LEFT); - x1 = findContourNeg (image, threshold, xR, yR, RIGHT); - } - // fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); - - xVec->data.F32[Npt + 0] = image->col0 + x0; - xVec->data.F32[Npt + 1] = image->col0 + x1; - yVec->data.F32[Npt + 0] = image->row0 + yR; - yVec->data.F32[Npt + 1] = image->row0 + yR; - Npt += 2; - - if (Npt >= xVec->nalloc - 1) { - psVectorRealloc (xVec, xVec->nalloc + 100); - psVectorRealloc (yVec, yVec->nalloc + 100); - } - yR --; - } -pt2: - xVec->n = Npt; - yVec->n = Npt; - - // fprintf (stderr, "done\n"); - psArray *tmpArray = psArrayAlloc(2); - tmpArray->data[0] = (psPtr *) xVec; - tmpArray->data[1] = (psPtr *) yVec; - psTrace(__func__, 3, "---- %s() end ----\n", __func__); - return(tmpArray); -} Index: unk/psphot/src/pmSourceFitSet.c =================================================================== --- /trunk/psphot/src/pmSourceFitSet.c (revision 6570) +++ (revision ) @@ -1,256 +1,0 @@ -# include "psphot.h" - -bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) { - - double x, y, Io; - - FILE *f = fopen (fitset, "r"); - if (f == NULL) return false; - - psArray *modelSet = psArrayAlloc (16); - modelSet->n = 0; - - while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) { - pmModel *model = pmModelAlloc (oneModel->type); - - for (psS32 i = 0; i < model->params->n; i++) { - model->params->data.F32[i] = oneModel->params->data.F32[i]; - model->dparams->data.F32[i] = oneModel->dparams->data.F32[i]; - } - model->params->data.F32[1] = Io; - model->params->data.F32[2] = x; - model->params->data.F32[3] = y; - psArrayAdd (modelSet, 16, model); - } - - pmSourceFitSet (source, modelSet, PSF); - - // write out positive object - psphotSaveImage (NULL, source->pixels, "object.fits"); - - // subtract object, leave local sky - for (int i = 0; i < modelSet->n; i++) { - pmModel *model = modelSet->data[i]; - pmSourceSubModel (source->pixels, source->mask, model, false, false); - - fprintf (stderr, "output parameters (obj %d):\n", i); - for (int n = 0; n < model->params->n; n++) { - fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]); - } - } - - // write out - psphotSaveImage (NULL, source->pixels, "resid.fits"); - psphotSaveImage (NULL, source->mask, "mask.fits"); - return true; -} - -# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 -# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 - -bool pmSourceFitSet (pmSource *source, - psArray *modelSet, - const bool PSF) -{ - psTrace(__func__, 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(source, false); - PS_ASSERT_PTR_NON_NULL(source->moments, false); - PS_ASSERT_PTR_NON_NULL(source->peak, false); - PS_ASSERT_PTR_NON_NULL(source->pixels, false); - PS_ASSERT_PTR_NON_NULL(source->mask, false); - PS_ASSERT_PTR_NON_NULL(source->weight, false); - - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - - // base values on first model - pmModel *model = modelSet->data[0]; - - // set the static variables - pmModelFitSetInit (model->type); - - int nSrc = modelSet->n; - int nPar = model->params->n - 1; // number of object parameters (excluding sky) - - // define parameter vectors for source set - psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - - pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); - - // define limits for single-source model - psVector *oneDelta; - psVector *oneParMin; - psVector *oneParMax; - modelLimits (&oneDelta, &oneParMin, &oneParMax); - - psMinConstrain *constrain = psMinConstrainAlloc(); - constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); - constrain->paramMask = PSF ? psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8) : NULL; - - // all but the sky are allowed to vary independently (subject to PSF) - // set the parameters from the multiple models - // also, set limits based on single-source limits - params->data.F32[0] = model->params->data.F32[0]; - constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0]; - constrain->paramMin->data.F32[0] = oneParMin->data.F32[0]; - constrain->paramMax->data.F32[0] = oneParMax->data.F32[0]; - constrain->paramMask->data.U8[0] = 0; - for (int i = 0; i < nSrc; i++) { - model = modelSet->data[i]; - for (int n = 1; n < nPar + 1; n++) { - params->data.F32[i*nPar + n] = model->params->data.F32[n]; - dparams->data.F32[i*nPar + n] = model->dparams->data.F32[n]; - constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n]; - constrain->paramMin->data.F32[i*nPar + n] = oneParMin->data.F32[n]; - constrain->paramMax->data.F32[i*nPar + n] = oneParMax->data.F32[n]; - if (PSF) { - constrain->paramMask->data.U8[i*nPar + n] = (n < 4) ? 0 : 1; - } - } - } - psFree (oneDelta); - psFree (oneParMin); - psFree (oneParMax); - - if (psTraceGetLevel(__func__) >= 5) { - for (int i = 0; i < params->n; i++) { - fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]); - } - } - - // PSF model only fits first 4 parameters, EXT model fits all - int nParams = PSF ? nSrc*3 + 1 : nSrc*nPar + 1; - - // maximum number of valid pixels - psS32 nPix = source->pixels->numRows * source->pixels->numCols; - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(nPix); - psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); - - nPix = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - // skip masked points - if (source->mask->data.U8[i][j]) { - continue; - } - // skip zero-weight points - if (source->weight->data.F32[i][j] == 0) { - continue; - } - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - // Convert i/j to image space: - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[nPix] = (psPtr *) coord; - y->data.F32[nPix] = source->pixels->data.F32[i][j]; - // psMinimizeLMChi2 takes wt = 1/dY^2 - yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; - nPix++; - } - } - if (nPix < nParams + 1) { - psTrace (__func__, 4, "insufficient valid pixels\n"); - psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__); - model->status = PM_MODEL_BADARGS; - psFree (x); - psFree (y); - psFree (yErr); - psFree (params); - psFree (dparams); - psFree(constrain->paramMask); - psFree(constrain->paramMin); - psFree(constrain->paramMax); - psFree(constrain->paramDelta); - psFree(constrain); - return(false); - } - x->n = nPix; - y->n = nPix; - yErr->n = nPix; - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - - psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); - - psTrace (__func__, 5, "fitting function\n"); - fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet); - - // parameter errors from the covariance matrix - for (int i = 0; i < dparams->n; i++) { - if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - - // get the Gauss-Newton distance for fixed model parameters - if (constrain->paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); - for (int i = 0; i < dparams->n; i++) { - if (!constrain->paramMask->data.U8[i]) - continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - psFree (delta); - } - - // assign back the parameters to the models - for (int i = 0; i < nSrc; i++) { - model = modelSet->data[i]; - model->params->data.F32[0] = params->data.F32[0]; - for (int n = 1; n < nPar + 1; n++) { - model->params->data.F32[n] = params->data.F32[i*nPar + n]; - model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; - } - // save the resulting chisq, nDOF, nIter - // these are not unique for any one source - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // set the model success or failure status - model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; - - // models can go insane: reject these - onPic &= (model->params->data.F32[2] >= source->pixels->col0); - onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (model->params->data.F32[3] >= source->pixels->row0); - onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) { - model->status = PM_MODEL_OFFIMAGE; - } - } - - source->mode |= PM_SOURCE_FITTED; - - psFree(x); - psFree(y); - psFree(yErr); - psFree(myMin); - psFree(covar); - psFree(constrain->paramMask); - psFree(constrain->paramMin); - psFree(constrain->paramMax); - psFree(constrain->paramDelta); - psFree(constrain); - psFree(params); - psFree(dparams); - - // free static memory used by pmModelFitSet - pmModelFitSetClear (); - - rc = (onPic && fitStatus); - psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); - psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc); - return(rc); -} - Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 6570) +++ /trunk/psphot/src/psModulesUtils.c (revision 6571) @@ -1,232 +1,3 @@ # include "psModulesUtils.h" - -// extract config informatin from config data or from image header, if specified -// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) { - - char *source; - char *keyword; - psF32 value; - psMetadataItem *item; - - // look for the entry in the config collection - item = psMetadataLookup (config, name); - if (item == NULL) { - psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name); - value = psMetadataLookupF32 (status, header, name); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in header", name); - } - *status = true; - return (value); - } - - // I'm either expecting a string, with the name "HD:keyword"... - if (item->type == PS_DATA_STRING) { - source = item->data.V; - if (!strncasecmp (source, "HD:", 3)) { - keyword = &source[3]; - value = psMetadataLookupF32 (status, header, keyword); - if (*status == false) { - psAbort ("pmConfigLookupF32", "no key %s in config", name); - } - *status = true; - // psFree (item); - return (value); - } - } - - // ... or a value (F32?) - if (item->type == PS_DATA_F32) { - value = item->data.F32; - // psFree (item); - return (value); - } - - *status = false; - psError(PS_ERR_UNKNOWN, true, "invalid item"); - return (0); -} - -bool pmModelFitStatus (pmModel *model) { - - bool status; - - pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type); - status = statusFunc (model); - - return (status); -} - -pmModel *pmModelSelect (pmSource *source) { - switch (source->type) { - case PM_SOURCE_STAR: - return source->modelPSF; - - case PM_SOURCE_EXTENDED: - return source->modelEXT; - - default: - return NULL; - } - return NULL; -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) { - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - *fitMag = 99.0; - *obsMag = 99.0; - - pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - if (fitSum <= 0) return false; - if (!isfinite(fitSum)) return false; - *fitMag = -2.5*log10(fitSum); - - for (int ix = 0; ix < image->numCols; ix++) { - for (int iy = 0; iy < image->numRows; iy++) { - if (mask->data.U8[iy][ix]) continue; - obsSum += image->data.F32[iy][ix] - sky; - } - } - if (obsSum <= 0) return false; - *obsMag = -2.5*log10(obsSum); - - return (true); -} - -// force the fitted sky to meet the source flux at outer radius -bool pmModelSkyOffset (pmModel *model, float x, float y, float radius) { - - float flux; - - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - - pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - psVector *params = model->params; - - coord->data.F32[0] = x + radius; - coord->data.F32[1] = y; - flux = modelFunc (NULL, params, coord); - params->data.F32[0] = flux; - - psFree (coord); - - return true; -} - -pmModel *pmModelCopy (pmModel *model) { - - pmModel *new = pmModelAlloc (model->type); - - new->chisq = model->chisq; - new->nDOF = model->nDOF; - new->nIter = model->nIter; - new->status = model->status; - new->radius = model->radius; - - for (int i = 0; i < new->params->n; i++) { - new->params->data.F32[i] = model->params->data.F32[i]; - new->dparams->data.F32[i] = model->dparams->data.F32[i]; - } - - return (new); -} - -float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) { - - int Xs, Xe, Ys, Ye; - int xi, xj, yi, yj; - int xIs, xJs, yIs, yJs; - int xIe, yIe; - float flux, wt; - - psImage *Pi = Mi->pixels; - psImage *Pj = Mj->pixels; - - psImage *Wi = Mi->weight; - - psImage *Ti = Mi->mask; - psImage *Tj = Mj->mask; - - Xs = PS_MAX (Pi->col0, Pj->col0); - Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); - - Ys = PS_MAX (Pi->row0, Pj->row0); - Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); - - xIs = Xs - Pi->col0; - xJs = Xs - Pj->col0; - yIs = Ys - Pi->row0; - yJs = Ys - Pj->row0; - - xIe = Xe - Pi->col0; - yIe = Ye - Pi->row0; - - // note that this is addressing the same image pixels, - // though only if both are source not model images - flux = 0; - for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { - for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.U8[yi][xi]) continue; - if (Tj->data.U8[yj][xj]) continue; - wt = Wi->data.F32[yi][xi]; - if (wt > 0) { - flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt; - } - } - } - return (flux); -} - -float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) { - - int Xs, Xe, Ys, Ye; - int xi, xj, yi, yj; - int xIs, xJs, yIs, yJs; - int xIe, yIe; - float flux, wt; - - psImage *Pi = Mi->pixels; - psImage *Pj = Mj->pixels; - - psImage *Wi = Mi->weight; - - psImage *Ti = Mi->mask; - psImage *Tj = Mj->mask; - - Xs = PS_MAX (Pi->col0, Pj->col0); - Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); - - Ys = PS_MAX (Pi->row0, Pj->row0); - Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); - - xIs = Xs - Pi->col0; - xJs = Xs - Pj->col0; - yIs = Ys - Pi->row0; - yJs = Ys - Pj->row0; - - xIe = Xe - Pi->col0; - yIe = Ye - Pi->row0; - - // note that this is addressing the same image pixels, - // though only if both are source not model images - flux = 0; - for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { - for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { - if (Ti->data.U8[yi][xi]) continue; - if (Tj->data.U8[yj][xj]) continue; - wt = Wi->data.F32[yi][xi]; - if (wt > 0) { - flux += 1.0 / wt; - } - } - } - return (flux); -} static void ppConfigFree (ppConfig *config) { Index: /trunk/psphot/src/psModulesUtils.h =================================================================== --- /trunk/psphot/src/psModulesUtils.h (revision 6570) +++ /trunk/psphot/src/psModulesUtils.h (revision 6571) @@ -23,5 +23,4 @@ psMetadata *recipe; // The recipe (i.e., specific setups) psMetadata *arguments; // Command-line arguments - psMetadata *options; // Command-line recipe options psDB *database; // Database handle } ppConfig; Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 6570) +++ /trunk/psphot/src/psphot.c (revision 6571) @@ -10,18 +10,14 @@ // load command-line arguments, options, and system config data - ppConfig *config = psphotArguments (&argc, argv); + pmConfig *config = psphotArguments (&argc, argv); // load input data (config and images (signal, noise, mask) - ppFile *input = psphotParseCamera (config); - - // check on output mode,format, setup basic filename, etc - psphotOutputPrep (input, config); + psphotParseCamera (config); // call psphot for each readout - psphotImageLoop (input, config); + psphotImageLoop (config); psLogMsg ("psphot", 3, "complete psphot run: %f sec\n", psTimerMark ("complete")); - psFree (input); psFree (config); psphotCleanup (); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 6570) +++ /trunk/psphot/src/psphot.h (revision 6571) @@ -10,7 +10,7 @@ # include # include -# include "psLibUtils.h" -# include "psModulesUtils.h" -# include "psSparse.h" +// # include "psLibUtils.h" +// # include "psModulesUtils.h" +// # include "psSparse.h" # include "pmFPAConstruct.h" # include "pmConcepts.h" @@ -21,12 +21,13 @@ // top-level psphot functions -ppConfig *psphotArguments (int *argc, char **argv); -ppFile *psphotParseCamera (ppConfig *config); -bool psphotImageLoop (ppFile *input, ppConfig *config); +pmConfig *psphotArguments (int *argc, char **argv); +bool psphotParseCamera (pmConfig *config); +bool psphotImageLoop (pmConfig *config); bool psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe); bool psphotReadout (pmReadout *readout, psMetadata *config); -bool ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum); void psphotCleanup (void); + +// bool ppImageLoadPixels (ppFile *input, psDB *db, int chipNum, int cellNum); // psphotReadout functions @@ -47,15 +48,5 @@ bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf); -// XXX deprecate these? -// bool psphotFullFit (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -// bool psphotApplyPSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -// bool psphotFitExtended (pmReadout *readout, psMetadata *config, psArray *sources, psStats *skyStats); - // basic support functions -pmModel *pmModelCopy (pmModel *model); -pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); -float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); -float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj); -psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); void psphotModelGroupInit (void); int psphotSortBySN (const void **a, const void **b); @@ -63,5 +54,5 @@ bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf); void psphotTestArguments (int *argc, char **argv); -bool pmCellSetMask (pmCell *cell, psMetadata *recipe); +bool psphotMaskCell (pmCell *cell, psMetadata *recipe); bool psphotBackgroundNames (psMetadata *arguments); bool psphotSkyReplace (pmReadout *readout, psImage *background); @@ -73,8 +64,7 @@ bool psphotInitRadiusEXT (psMetadata *config, pmModelType type); bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model); -bool psphotDefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); -bool psphotRedefinePixels (pmSource *mySource, const pmReadout *readout, psF32 x, psF32 y, psF32 Radius); // output functions +# if (0) bool pmSourcesWriteSX (psArray *sources, char *filename); bool pmSourcesWriteOBJ (psArray *sources, char *filename); @@ -93,11 +83,8 @@ int psphotSaveImage (psMetadata *header, psImage *image, char *filename); bool psphotUpdateHeader (psMetadata *header, psMetadata *config); -int pmSourcesDophotType (pmSource *source); -bool psMetadataItemTransfer (psMetadata *out, psMetadata *in, char *key); +# endif -char *psphotSplitName (psMetadata *header); -void psphotOutputPrep (ppFile *file, ppConfig *config); -void psphotOutputCleanup (void); -char *psphotNameSubstitute (char *input, char *replace, char *key); +// void psphotOutputPrep (ppFile *file, ppConfig *config); +// void psphotOutputCleanup (void); bool psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf); @@ -107,11 +94,4 @@ bool psphotEvalDBL (pmSource *source, pmModel *model); bool psphotEvalEXT (pmSource *source, pmModel *model); - -// functions to support simultaneous multi-source fitting -bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); -bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); -psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); -bool pmModelFitSetInit (pmModelType type); -void pmModelFitSetClear (void); // functions to support the source fitting process @@ -124,10 +104,5 @@ psArray *psphotFitDBL (pmReadout *readout, pmSource *source); -// bicubic interpolation -psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); -psPlane psImageBicubeMin (psPolynomial2D *poly); - -bool psImageJpegColormap (char *name); -bool psImageJpeg (psImage *image, char *filename, float zero, float scale); +// XXX these can probably be dropped: // optional mode for clip fit? @@ -143,2 +118,26 @@ const psVector *z, const psVector *t); + +// functions to support simultaneous multi-source fitting +bool psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF); +bool pmSourceFitSet (pmSource *source, psArray *modelSet, const bool PSF); +psF32 pmModelFitSet (psVector *deriv, const psVector *params, const psVector *x); +bool pmModelFitSetInit (pmModelType type); +void pmModelFitSetClear (void); + +// bicubic interpolation +psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y); +psPlane psImageBicubeMin (psPolynomial2D *poly); + +bool psImageJpegColormap (char *name); +bool psImageJpeg (psImage *image, char *filename, float zero, float scale); + + +// XXX deprecate + +// char *psphotNameSubstitute (char *input, char *replace, char *key); +// char *psphotSplitName (psMetadata *header); +// pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius); +// float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj); +// float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj); +// psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 6570) +++ /trunk/psphot/src/psphotArguments.c (revision 6571) @@ -1,7 +1,8 @@ # include "psphot.h" +# include static void usage (void); -ppConfig *psphotArguments (int *argc, char **argv) { +pmConfig *psphotArguments (int *argc, char **argv) { int N; @@ -12,50 +13,97 @@ psLogSetFormat ("M"); - ppConfig *config = ppConfigAlloc (); + // these other options override the PSPHOT recipe options + psMetadata *options = psMetadataAlloc (); - // load config data from default locations - if (!pmConfigRead(&config->site, &config->camera, &config->recipe, argc, argv, PSPHOT_RECIPE)) { - psErrorStackPrint(stderr, "Can't find site configuration!\n"); - exit (EXIT_FAILURE); + // run the test model (requires X,Y coordinate) + if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { + psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT", 0, "", true); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); + psMetadataAddF32 (options, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); + + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + + // specify the modeltest model + if ((N = psArgumentGet (*argc, argv, "-model"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]); + psArgumentRemove (N, argc, argv); + } + + // specify the test fit mode + if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]); + psArgumentRemove (N, argc, argv); + } + if ((N = psArgumentGet (*argc, argv, "-fitset"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]); + psArgumentRemove (N, argc, argv); + } } - // Parse other command-line arguments - config->arguments = psMetadataAlloc (); - - // arguments (must be supplied for each run, or not used) - // mask image (input) is used by psphotImageLoop - if ((N = psArgumentGet (*argc, argv, "-mask"))) { + // photcode : used in output to supplement header data (argument or recipe?) + if ((N = psArgumentGet (*argc, argv, "-photcode"))) { psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "MASK_IMAGE", PS_META_REPLACE, "", argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, argc, argv); } - // weight image (input) is used by psphotImageLoop - if ((N = psArgumentGet (*argc, argv, "-weight"))) { + // break : used from recipe throughout psphotReadout + if ((N = psArgumentGet (*argc, argv, "-break"))) { psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "WEIGHT_IMAGE", PS_META_REPLACE, "", argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, argc, argv); } - // residual image (output) is used by psphotOutput - if ((N = psArgumentGet (*argc, argv, "-resid"))) { + // fitmode : used from recipe throughout psphotReadout + if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "RESID_IMAGE", PS_META_REPLACE, "", argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, argc, argv); } - // background image (output) is used by psphotOutput - if ((N = psArgumentGet (*argc, argv, "-background"))) { + // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF + if ((N = psArgumentGet (*argc, argv, "-region"))) { psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_IMAGE", PS_META_REPLACE, "", argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); psArgumentRemove (N, argc, argv); } - // background model (output) is used by psphotOutput - if ((N = psArgumentGet (*argc, argv, "-backmodel"))) { + // other arbitrary recipe values: -D key value (all added as string) + while ((N = psArgumentGet (*argc, argv, "-D"))) { psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "BACKGROUND_MODEL", PS_META_REPLACE, "", argv[N]); + psMetadataAddStr (options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); + psArgumentRemove (N, argc, argv); psArgumentRemove (N, argc, argv); } + + // other arbitrary recipe values: -Df key value (all added as float) + while ((N = psArgumentGet (*argc, argv, "-Df"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddF32 (options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // other arbitrary recipe values: -Di key value (all added as int) + while ((N = psArgumentGet (*argc, argv, "-Di"))) { + psArgumentRemove (N, argc, argv); + psMetadataAddS32 (options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + + // load config data from default locations + pmConfig *config = pmConfigRead(argc, argv, PSPHOT_RECIPE); + + // Storage for other command-line arguments + config->arguments = psMetadataAlloc (); + + // save these recipe options until we have loaded the options + psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "PSPHOT.OPTIONS", PS_META_REPLACE, "", options); // chip selection is used to limit chips to be processed @@ -66,108 +114,58 @@ } - // run the test model (requires X,Y coordinate) - if ((N = psArgumentGet (*argc, argv, "-modeltest"))) { - psMetadataAddBool (config->arguments, PS_LIST_TAIL, "TEST_FIT", 0, "", true); - psMetadataAddF32 (config->arguments, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1])); - psMetadataAddF32 (config->arguments, PS_LIST_TAIL, "TEST_FIT_Y", 0, "", atof(argv[N+2])); + // we load all input files onto a psArray, to be parsed later + psArray *input = psArrayAlloc (16); + // load the list of filenames the supplied file (may be a glob: "file*.fits") + if ((N = psArgumentGet (*argc, argv, "-file"))) { + glob_t globList; psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - - // specify the modeltest model - if ((N = psArgumentGet (*argc, argv, "-model"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]); - psArgumentRemove (N, argc, argv); + globList.gl_offs = 0; + glob (argv[N], 0, NULL, &globList); + for (int i = 0; i < globList.gl_pathc; i++) { + char *filename = psStringCopy (globList.gl_pathv[i]); + psArrayAdd (input, 16, filename); } - - // specify the test fit mode - if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]); - psArgumentRemove (N, argc, argv); - } - if ((N = psArgumentGet (*argc, argv, "-fitset"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]); - psArgumentRemove (N, argc, argv); - } - } - - // these other options override the recipe options - config->options = psMetadataAlloc (); - - // PSF : optional psf file to be loaded - if ((N = psArgumentGet (*argc, argv, "-psf"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_INPUT_FILE", PS_META_REPLACE, "", argv[N]); psArgumentRemove (N, argc, argv); } - // PSF : optional psf file to be loaded - if ((N = psArgumentGet (*argc, argv, "-psfout"))) { + // load the list from the supplied text file + if ((N = psArgumentGet (*argc, argv, "-list"))) { + int nItems; + char line[1024]; // XXX limits the list lines to 1024 chars + char word[1024]; + char *filename; + psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "PSF_OUTPUT_FILE", PS_META_REPLACE, "", argv[N]); + FILE *f = fopen (argv[N], "r"); + if (f == NULL) { + psAbort ("psphot", "unable to open specified list file"); + } + while (fgets (line, 1024, f) != NULL) { + nItems = sscanf (line, "%s", word); + switch (nItems) { + case 0: + break; + case 1: + filename = psStringCopy (word); + psArrayAdd (input, 16, filename); + break; + default: + // rigid format, no comments allowed? + psAbort ("psphot", "error parsing input list file"); + break; + } + } psArgumentRemove (N, argc, argv); - } + } + if (input->n < 1) usage (); - // photcode : used in output to supplement header data (argument or recipe?) - if ((N = psArgumentGet (*argc, argv, "-photcode"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); - psArgumentRemove (N, argc, argv); - } + // input list gets places as an array on the config->arguements list + psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "INPUT", PS_META_REPLACE, "", input); - // break : used from recipe throughout psphotReadout - if ((N = psArgumentGet (*argc, argv, "-break"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]); - psArgumentRemove (N, argc, argv); - } - - // fitmode : used from recipe throughout psphotReadout - if ((N = psArgumentGet (*argc, argv, "-fitmode"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]); - psArgumentRemove (N, argc, argv); - } - - // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF - if ((N = psArgumentGet (*argc, argv, "-region"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]); - psArgumentRemove (N, argc, argv); - } - - // other arbitrary recipe values: -D key value (all added as string) - while ((N = psArgumentGet (*argc, argv, "-D"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddStr (config->options, PS_LIST_TAIL, argv[N], 0, "", argv[N+1]); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - - // other arbitrary recipe values: -Df key value (all added as float) - while ((N = psArgumentGet (*argc, argv, "-Df"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddF32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atof(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - - // other arbitrary recipe values: -Di key value (all added as int) - while ((N = psArgumentGet (*argc, argv, "-Di"))) { - psArgumentRemove (N, argc, argv); - psMetadataAddS32 (config->options, PS_LIST_TAIL, argv[N], 0, "", atoi(argv[N+1])); - psArgumentRemove (N, argc, argv); - psArgumentRemove (N, argc, argv); - } - - if (*argc != 3) usage (); + if (*argc != 2) usage (); // input and output positions are fixed - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "INPUT_FILE", PS_META_REPLACE, "", argv[1]); - psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT_ROOT", PS_META_REPLACE, "", argv[2]); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", PS_META_REPLACE, "", argv[1]); psTrace(__func__, 1, "Done with psphotArguments...\n"); @@ -176,5 +174,5 @@ static void usage (void) { - fprintf (stderr, "USAGE: psphot (image) (output)\n"); + fprintf (stderr, "USAGE: psphot [-file image(s)] [-list imagelist] (output)\n"); exit (2); } Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 6570) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 6571) @@ -3,5 +3,5 @@ // 2006.02.07 : no leaks! // try PSF models and select best option -pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources) { +pmPSF *psphotChoosePSF (psArray *sources, psMetadata *config) { bool status; Index: /trunk/psphot/src/psphotCleanup.c =================================================================== --- /trunk/psphot/src/psphotCleanup.c (revision 6570) +++ /trunk/psphot/src/psphotCleanup.c (revision 6571) @@ -9,5 +9,5 @@ psTimeFinalize (); pmConceptsDone (); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot"); + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot"); return; } Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 6570) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 6571) @@ -4,5 +4,5 @@ // 2006.02.07 : no leaks! // fit all reasonable sources with the linear PSF model -bool psphotEnsemblePSF (pmReadout *readout, psMetadata *config, psArray *sources, pmPSF *psf, bool final) { +bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf, bool final) { bool status; Index: /trunk/psphot/src/psphotFitSet.c =================================================================== --- /trunk/psphot/src/psphotFitSet.c (revision 6571) +++ /trunk/psphot/src/psphotFitSet.c (revision 6571) @@ -0,0 +1,47 @@ +# include "psphot.h" + +bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) { + + double x, y, Io; + + FILE *f = fopen (fitset, "r"); + if (f == NULL) return false; + + psArray *modelSet = psArrayAlloc (16); + modelSet->n = 0; + + while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) { + pmModel *model = pmModelAlloc (oneModel->type); + + for (psS32 i = 0; i < model->params->n; i++) { + model->params->data.F32[i] = oneModel->params->data.F32[i]; + model->dparams->data.F32[i] = oneModel->dparams->data.F32[i]; + } + model->params->data.F32[1] = Io; + model->params->data.F32[2] = x; + model->params->data.F32[3] = y; + psArrayAdd (modelSet, 16, model); + } + + pmSourceFitSet (source, modelSet, PSF); + + // write out positive object + psphotSaveImage (NULL, source->pixels, "object.fits"); + + // subtract object, leave local sky + for (int i = 0; i < modelSet->n; i++) { + pmModel *model = modelSet->data[i]; + pmSourceSubModel (source->pixels, source->mask, model, false, false); + + fprintf (stderr, "output parameters (obj %d):\n", i); + for (int n = 0; n < model->params->n; n++) { + fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]); + } + } + + // write out + psphotSaveImage (NULL, source->pixels, "resid.fits"); + psphotSaveImage (NULL, source->mask, "mask.fits"); + return true; +} + Index: /trunk/psphot/src/psphotImageMedian.c =================================================================== --- /trunk/psphot/src/psphotImageMedian.c (revision 6570) +++ /trunk/psphot/src/psphotImageMedian.c (revision 6571) @@ -9,32 +9,25 @@ static int MAX_SAMPLE_PIXELS; -static char *backImage = NULL; -static char *backModel = NULL; - -bool psphotBackgroundNames (psMetadata *arguments) { - - bool status; - - backImage = psMetadataLookupStr (&status, arguments, "BACKGROUND_IMAGE"); - backModel = psMetadataLookupStr (&status, arguments, "BACKGROUND_MODEL"); - - return true; -} - // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -psImage *psphotImageMedian (pmReadout *readout, psMetadata *config) +psImage *psphotImageMedian (psMetadata *config, pmFPAview *view) { bool status; psRegion region; + psImage *model = NULL; + + psTimerStart ("psphot"); + + // find the currently selected readout + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); + pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + pmReadout *readout = pmFPAviewThisReadout (view, input); + + psImage *image = readout->image; + psImage *mask = readout->mask; rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, config, "IMSTATS_NPIX"); if (!status) MAX_SAMPLE_PIXELS = 1000; - - psTimerStart ("psphot"); - - psImage *image = readout->image; - psImage *mask = readout->mask; // dimensions of input & output image @@ -56,6 +49,6 @@ int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY; - // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - psImage *model = psImageAlloc (nx, ny, PS_TYPE_F32); + // select model pixels, from output background model file, or create + model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32); // measure clipped median for subimages @@ -76,6 +69,7 @@ // XXX psImageStats is still very inefficient and poorly coded... + // psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); // stats = psImageStats (stats, subset, maskset, 0xff); - // model->data.F32[iy][ix] = stats->clippedMean; + // backMdl->data.F32[iy][ix] = stats->clippedMean; } } @@ -83,6 +77,7 @@ // linear interpolation to full-scale - - psImage *background = psImageAlloc (Nx, Ny, PS_TYPE_F32); + + // select background pixels, from output background file, or create + background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32); // XXX this code skips the initial pixels @@ -216,4 +211,7 @@ psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot")); + // back-sub image pixels, from output background file (don't create if not requested) + backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32); + // subtract the background model for (int j = 0; j < image->numRows; j++) { @@ -221,16 +219,19 @@ if (!mask->data.U8[j][i]) { image->data.F32[j][i] -= background->data.F32[j][i]; - } - } - } - - // optionally save background - if (backImage != NULL) psphotSaveImage (NULL, background, backImage); - if (backModel != NULL) psphotSaveImage (NULL, model, backModel); + if (backSub) { + backSub->data.F32[j][i] = image->data.F32[j][i]; + } + } + } + } psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot")); psFree (rnd); + + // it is safe to free all of these: either they are saved on config->files or are local psFree (model); - return (background); + psFree (background); + psFree (backSub); + return true; } Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 6570) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 6571) @@ -1,74 +1,3 @@ # include "psphot.h" - -// XXX EAM : the apMag should only be calculated for the brighter sources? -// XXX EAM : SN limit set by user? -// XXX EAM : masked region should be (optionally) elliptical -pmModel *pmSourceMagnitudes (pmSource *source, pmPSF *psf, float apRadius) { - - int status; - bool isPSF; - float x, y; - float rflux; - float radius; - pmModel *model; - - switch (source->type) { - case PM_SOURCE_STAR: - model = source->modelPSF; - if (model == NULL) return NULL; - radius = (apRadius > 0) ? apRadius : model->radius; - isPSF = true; - break; - - case PM_SOURCE_EXTENDED: - model = source->modelEXT; - if (model == NULL) return NULL; - radius = model->radius; - isPSF = false; - break; - - default: - return NULL; - } - - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - - // replace source flux - pmSourceAddModel (source->pixels, source->mask, model, false, false); - - // set aperture mask circle of PSF_FIT_RADIUS - psImageKeepCircle (source->mask, x, y, radius, "OR", PSPHOT_MASK_MARKED); - - // measure object photometry - status = pmSourcePhotometry (&source->fitMag, &source->apMag, model, source->pixels, source->mask); - - // for PSFs, correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS - if (isPSF && (psf != NULL)) { - if (psf->growth != NULL) { - source->apMag += pmGrowthCurveCorrect (psf->growth, model->radius); - } - - rflux = pow (10.0, 0.4*source->fitMag); - source->apMag -= PS_SQR(model->radius)*rflux * psf->skyBias + psf->skySat / rflux; - source->fitMag += psPolynomial4DEval (psf->ApTrend, x, y, 0.0, 0.0); - } - - // unmask aperture - psImageKeepCircle (source->mask, x, y, radius, "AND", ~PSPHOT_MASK_MARKED); - - // subtract object, leave local sky - pmSourceSubModel (source->pixels, source->mask, model, false, false); - - if (!status) return NULL; - return model; -} - -/* -aprMag' - fitMag = flux*skySat + r^2*rflux*skyBias + ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) - fitMAg = ApTrend(x,y) -(aprMag - flux*skySat - r^2*rflux*skyBias) = fitMAg + ApTrend(x,y) - -*/ bool psphotMagnitudes (psMetadata *config, psArray *sources, pmPSF *psf) { @@ -88,2 +17,3 @@ } +// XXX should this function use RADIUS or should it Index: /trunk/psphot/src/psphotMaskCell.c =================================================================== --- /trunk/psphot/src/psphotMaskCell.c (revision 6571) +++ /trunk/psphot/src/psphotMaskCell.c (revision 6571) @@ -0,0 +1,31 @@ +# include "psphot.h" + +// 2006.02.04 : no leaks +bool psphotMaskCell (pmCell *cell, psMetadata *recipe) { + + bool status; + + // mask the excluded outer pixels + // these coordinates refer to the parent image + // these bounds will saturate on the subimage + // negative upper bounds will subtract from the *subimage* + float XMIN = psMetadataLookupF32 (&status, recipe, "XMIN"); + float XMAX = psMetadataLookupF32 (&status, recipe, "XMAX"); + float YMIN = psMetadataLookupF32 (&status, recipe, "YMIN"); + float YMAX = psMetadataLookupF32 (&status, recipe, "YMAX"); + psRegion valid = psRegionSet (XMIN, XMAX, YMIN, YMAX); + + // Set the pixels + psArray *readouts = cell->readouts; // Array of readouts + for (int i = 0; i < readouts->n; i++) { + pmReadout *readout = readouts->data[i]; // The readout of interest + + // psRegionForImage creates a region for the subimage, not the parent + // XXX I propose that psRegionForImage treat the region as referring to the parent... + psRegion keep = psRegionForImage (readout->image, valid); + + // psImageKeepRegion assumes the region refers to the parent coordinates + psImageKeepRegion (readout->mask, keep, "OR", PSPHOT_MASK_INVALID); + } + return true; +} Index: /trunk/psphot/src/psphotModelGroupInit.c =================================================================== --- /trunk/psphot/src/psphotModelGroupInit.c (revision 6570) +++ /trunk/psphot/src/psphotModelGroupInit.c (revision 6571) @@ -1,3 +1,6 @@ # include "psphot.h" + +// Add locally-defined models here. As these mature, they can be moved to +// psModule/src/objects/models # include "models/pmModel_TAUSS.c" @@ -6,5 +9,4 @@ {"PS_MODEL_TAUSS", 9, pmModelFunc_TAUSS, pmModelFlux_TAUSS, pmModelRadius_TAUSS, pmModelLimits_TAUSS, pmModelGuess_TAUSS, pmModelFromPSF_TAUSS, pmModelFitStatus_TAUSS}, }; - void psphotModelGroupInit (void) Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 6570) +++ /trunk/psphot/src/psphotModelTest.c (revision 6571) @@ -3,5 +3,5 @@ static char DEFAULT_MODE[] = "EXT"; -bool psphotModelTest (pmReadout *readout, psMetadata *arguments, psMetadata *recipe) { +bool psphotModelTest (pmReadout *readout, recipe) { bool status; @@ -15,8 +15,8 @@ // run model fitting tests on a single source - if (!psMetadataLookupBool (&status, arguments, "TEST_FIT")) return false; + if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false; // what fitting mode to use? - char *psfModeWord = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODE"); + char *psfModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE"); if (!status) { psfModeWord = DEFAULT_MODE; @@ -33,5 +33,5 @@ } else { // find the model: supplied by user or first in the PSF_MODEL list - char *modelName = psMetadataLookupStr (&status, arguments, "TEST_FIT_MODEL"); + char *modelName = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODEL"); if (modelName == NULL) { // get the list pointers for the PSF_MODEL entries @@ -57,20 +57,20 @@ // find the fitting parameters (try test values first) - float INNER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_INNER_RADIUS"); + float INNER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_INNER_RADIUS"); if (!status) { INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); } - float OUTER = psMetadataLookupF32 (&status, arguments, "TEST_FIT_OUTER_RADIUS"); + float OUTER = psMetadataLookupF32 (&status, recipe, "TEST_FIT_OUTER_RADIUS"); if (!status) { OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); } - float RADIUS = psMetadataLookupF32 (&status, arguments, "TEST_FIT_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, recipe, "TEST_FIT_RADIUS"); if (!status) { RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); } - float mRADIUS = psMetadataLookupF32 (&status, arguments, "TEST_MOMENTS_RADIUS"); + float mRADIUS = psMetadataLookupF32 (&status, recipe, "TEST_MOMENTS_RADIUS"); if (!status) { mRADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); @@ -78,6 +78,6 @@ // define the source of interest - float xObj = psMetadataLookupF32 (&status, arguments, "TEST_FIT_X"); - float yObj = psMetadataLookupF32 (&status, arguments, "TEST_FIT_Y"); + float xObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X"); + float yObj = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y"); // construct the source structures @@ -120,5 +120,5 @@ } sprintf (name, "TEST_FIT_PAR%d", i); - value = psMetadataLookupF32 (&status, arguments, name); + value = psMetadataLookupF32 (&status, recipe, name); if (status) { params[i] = value; @@ -154,5 +154,5 @@ psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED); - char *fitset = psMetadataLookupStr (&status, arguments, "TEST_FIT_SET"); + char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); if (status) { status = psphotFitSet (source, model, fitset, psfMode); Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 6570) +++ /trunk/psphot/src/psphotOutput.c (revision 6571) @@ -1,153 +1,3 @@ # include "psphot.h" - -static int extNumber = -1; - -static char *outputName = NULL; -static char *outputRoot = NULL; -static char *outputMode = NULL; -static char *outputFormat = NULL; - -static char *extNumberFormat = NULL; -static char *extNameKey = NULL; -static char *extRoot = NULL; -static char *psfFile = NULL; -static char *psfSample = NULL; - -static psFits *outputFits = NULL; -static FILE *outputFile = NULL; - -bool psphotDataIO (psphotData *data, ppConfig *config, int chip, int cell) { - - // load input data, if depth is appropriate - ppImageLoadDepth loadDepth = ppImageCheckDepth (config->recipe, "LOAD.DEPTH"); - if ((chip == -1) && (cell == -1) && (loadDepth == PP_LOAD_FPA)) { - psTrace(__func__, 1, "Loading pixels for FPA...\n"); - ppImageLoadPixels(data->input, config->database, -1, -1); - } - if ((chip != -1) && (cell == -1) && (loadDepth == PP_LOAD_CHIP)) { - psTrace(__func__, 1, "Loading pixels for chip %d...\n", i); - ppImageLoadPixels(data->input, config->database, chip, -1); - } - if ((chip != -1) && (cell != -1) && (loadDepth == PP_LOAD_CELL)) { - psTrace(__func__, 1, "Loading pixels for chip %d, cell %d...\n", i, j); - ppImageLoadPixels(data->input, config->database, chip, cell); - } - - - - if (saveDepth == PP_LOAD_FPA) { - psTrace(__func__, 1, "Saving objects for FPA...\n"); - psphotOutputOpen (file->fpa->phu); - } - - ppImageLoadDepth saveDepth = ppImageCheckDepth (config->recipe, "SAVE.DEPTH"); - - -} - -void psphotOutputPrep (ppFile *file, ppConfig *config) { - - bool status; - - outputRoot = psMetadataLookupStr (&status, config->arguments, "OUTPUT_ROOT"); - if (!status) psAbort ("psphot", "output file not specified"); - - outputName = psMetadataLookupStr (&status, config->recipe, "OUTPUT_NAME"); - outputMode = psMetadataLookupStr (&status, config->recipe, "OUTPUT_MODE"); - outputFormat = psMetadataLookupStr (&status, config->recipe, "OUTPUT_FORMAT"); - - psfFile = psMetadataLookupStr (&status, config->recipe, "PSF_OUTPUT_FILE"); - psfSample = psMetadataLookupStr (&status, config->recipe, "PSF_SAMPLE_FILE"); - - // rules to construct output names - extNumberFormat = psMetadataLookupStr (&status, config->recipe, "EXTNUMBER_FORMAT"); - extNameKey = psMetadataLookupStr (&status, config->recipe, "EXTNAME"); - extRoot = psMetadataLookupStr (&status, config->recipe, "EXTROOT"); - - headExtname = psMetadataLookupStr (&status, config->recipe, "HEAD_EXTNAME"); - dataExtname = psMetadataLookupStr (&status, config->recipe, "DATA_EXTNAME"); - - if (extNumberFormat == NULL) { - extNumberFormat = psStringCopy ("%02d"); - } - - psStringStrip (outputName); - psStringStrip (extNameKey); - psStringStrip (extRoot); - - // validate the outputMode and outputFormat - status = false; - status |= !strcasecmp (outputFormat, "TEXT"); - status |= !strcasecmp (outputFormat, "OBJ"); - status |= !strcasecmp (outputFormat, "SX"); - status |= !strcasecmp (outputFormat, "CMP"); - status |= !strcasecmp (outputFormat, "CMF"); - if (!status) { - psLogMsg ("psphot", PS_LOG_ERROR, "invalid output format %s", outputFormat); - exit (1); - } - - // validate the outputMode and outputFormat - status = false; - status |= !strcasecmp (outputMode, "MEF"); - status |= !strcasecmp (outputMode, "SPLIT"); - if (!status) { - psLogMsg ("psphot", PS_LOG_ERROR, "invalid output mode %s", outputMode); - exit (1); - } - - // for MEF output, we only allow certain output formats - if (!strcasecmp (outputMode, "MEF")) { - if (!strcasecmp (outputFormat, "CMF")) goto valid_format; - psLogMsg ("psphot", PS_LOG_ERROR, "MEF output mode unavailable"); - exit (1); - } -valid_format: - - // register background image-file names - psphotBackgroundNames (config->arguments); - - // save so freeing config will not drop these references - psMemCopy (outputRoot); - psMemCopy (outputName); - psMemCopy (outputMode); - psMemCopy (outputFormat); - psMemCopy (extNumberFormat); - psMemCopy (extNameKey); - psMemCopy (extRoot); - psMemCopy (psfFile); - psMemCopy (psfSample); - - return; -} - -bool psphotOutputOpen (psMetadata *phu){ - - /* outputFile is used by text-format outputs */ - - psFree (outputFile); - outputFile = psphotOutputName (phu, outputName); - - /* for fits-format output, open outputFits */ - if (!strcasecmp (outputFormat, "CMF")) { - if (outputFits != NULL) { - psFitsClose (outputFits); - } - outputFits = psFitsOpen (outputFile, "w"); - - if (!strcasecmp (outputMode, "MEF")) { - psFitsWriteHeader (outputFits, phu); - } - } - return true; -} - -bool psphotOutputClose (void) { - - if (outputFile != NULL) { - psFitsClose (outputFile); - } - return true; -} // output functions: we have several fixed modes (sx, obj, cmp) @@ -159,26 +9,8 @@ psMetadata *header = pmReadoutGetHeader (readout); - psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES"); pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - - // starting value is -1 - extNumber ++; - - // for SPLIT output, we need to open a file for each output call - if (!strcasecmp (outputMode, "SPLIT")) { - outputFile = psphotSplitName (header); - } - - if (!strcasecmp (outputMode, "MEF")) { - headExt = psphotOutputName (header, headExtname); - dataExt = psphotOutputName (header, dataExtname); - } - if (!strcasecmp (outputMode, "SPLIT")) { - dataExt = psphotOutputName (header, dataExtname); - } - - char *residImage = psMetadataLookupStr (&status, arguments, "RESID_IMAGE"); + // sample PSF images?? if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); - if (residImage != NULL) psphotSaveImage (header, readout->image, residImage); + if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); if (psfFile != NULL) { @@ -187,567 +19,4 @@ psFree (psfData); } - if (outputFile == NULL) { - psLogMsg ("output", 3, "no data output file selected"); - return; - } - if (outputFormat == NULL) { - psLogMsg ("output", 3, "no data output format selected"); - psFree (outputFile); - return; - } - if (!strcasecmp (outputFormat, "SX")) { - pmSourcesWriteSX (sources, outputFile); - psFree (outputFile); - return; - } - if (!strcasecmp (outputFormat, "OBJ")) { - pmSourcesWriteOBJ (sources, outputFile); - psFree (outputFile); - return; - } - if (!strcasecmp (outputFormat, "CMP")) { - pmSourcesWriteCMP (sources, outputFile, header); - psFree (outputFile); - return; - } - if (!strcasecmp (outputFormat, "CMF")) { - pmSourcesWriteCMF (sources, header, headExt, dataExt); - psFree (outputFile); - return; - } - if (!strcasecmp (outputFormat, "TEXT")) { - pmSourcesWriteText (sources, outputFile); - psFree (outputFile); - return; - } - - psAbort ("psphot", "unknown output mode %s", outputFormat); -} - -// dophot-style output list with fixed line width -bool pmSourcesWriteOBJ (psArray *sources, char *filename) { - - int type; - psF32 *PAR, *dPAR; - float dmag, apResid; - - psTimerStart ("string"); - - psLine *line = psLineAlloc (104); // 104 is dophot-defined line length - - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - - // write sources with models - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - dmag = dPAR[1] / PAR[1]; - type = pmSourceDophotType (source); - apResid = source->apMag - source->fitMag; - - psLineInit (line); - psLineAdd (line, "%3d", type); - psLineAdd (line, "%8.2f", PAR[2]); - psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", source->fitMag); - psLineAdd (line, "%6.3f", dmag); - psLineAdd (line, "%9.2f", PAR[0]); - psLineAdd (line, "%9.3f", PAR[4]); - psLineAdd (line, "%9.3f", PAR[5]); - psLineAdd (line, "%7.2f", PAR[6]); - psLineAdd (line, "%8.3f", 99.999); - psLineAdd (line, "%8.3f", source->apMag); - psLineAdd (line, "%8.2f\n", apResid); - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - psFree (line); - fprintf (stderr, "%f seconds for %d objects with psLine\n", psTimerMark ("string"), (int)sources->n); - - psTimerStart ("string"); - - f = fopen ("test.obj", "w"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", "test.obj"); - return false; - } - - char *string; - // write sources with models - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - dmag = dPAR[1] / PAR[1]; - type = pmSourceDophotType (source); - apResid = source->apMag - source->fitMag; - - string = NULL; - psStringAppend (&string, "%3d", type); - psStringAppend (&string, "%8.2f", PAR[2]); - psStringAppend (&string, "%8.2f", PAR[3]); - psStringAppend (&string, "%8.3f", source->fitMag); - psStringAppend (&string, "%6.3f", dmag); - psStringAppend (&string, "%9.2f", PAR[0]); - psStringAppend (&string, "%9.3f", PAR[4]); - psStringAppend (&string, "%9.3f", PAR[5]); - psStringAppend (&string, "%7.2f", PAR[6]); - psStringAppend (&string, "%8.3f", 99.999); - psStringAppend (&string, "%8.3f", source->apMag); - psStringAppend (&string, "%8.2f\n", apResid); - fwrite (string, 1, strlen(string), f); - psFree (string); - } - fclose (f); - fprintf (stderr, "%f seconds for %d objects with psString\n", psTimerMark ("string"), (int)sources->n); - - return true; -} - -// elixir-mode / sextractor-style output list with fixed line width -bool pmSourcesWriteSX (psArray *sources, char *filename) { - - psF32 *PAR, *dPAR; - float dmag; - - psLine *line = psLineAlloc (110); // 110 is sextractor line length - - FILE *f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("WriteSourceSX", 3, "can't open output file for output %s\n", filename); - return false; - } - - // write sources with models - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - // pmSourceSextractType (source, &type, &flags); XXX EAM : implement me... - - dmag = dPAR[1] / PAR[1]; - - psLineInit (line); - psLineAdd (line, "%5.2f", 0.0); // should be type - psLineAdd (line, "%11.3f", PAR[2]); - psLineAdd (line, "%11.3f", PAR[3]); - psLineAdd (line, "%9.4f", source->fitMag); - psLineAdd (line, "%9.4f", dmag); - psLineAdd (line, "%13.4f", PAR[0]); - psLineAdd (line, "%9.2f", 0.0); // should be FWHMx - psLineAdd (line, "%9.2f", 0.0); // should be FWHMy - psLineAdd (line, "%6.1f", 0.0); // should be Theta - psLineAdd (line, "%9.4f", 99.999); // should be MAG_ISO - psLineAdd (line, "%9.4f", source->apMag); - psLineAdd (line, "%4d\n", 0); // should be flags - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - return true; -} - -// elixir-style pseudo FITS table (header + ascii list) -bool pmSourcesWriteCMP (psArray *sources, char *filename, psMetadata *header) { - - int i, type; - psMetadataItem *mdi; - psF32 *PAR, *dPAR; - float dmag, lsky; - bool status; - - // find config information for output header - float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); - if (!status) ZERO_POINT = 25.0; - - // write necessary information to output header - psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", sources->n); - - // create file, write-out header - psFits *fits = psFitsOpen (filename, "w"); - - // set NAXIS to 0 : CFITSIO requires isolated header to have NAXIS = 0 - mdi = psMetadataLookup (header, "NAXIS"); - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; - - // psMetadataAdd (header, PS_LIST_HEAD, "NAXIS", PS_DATA_S32 | PS_META_REPLACE, "", 0); - psFitsWriteHeader (header, fits); - psFitsClose (fits); - - // re-open, add data to end of file - FILE *f = fopen (filename, "a+"); - if (f == NULL) { - psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename); - return false; - } - fseek (f, 0, SEEK_END); - - psLine *line = psLineAlloc (67); // 66 is imclean-defined line length - - // write sources with models first - for (i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - dmag = dPAR[1] / PAR[1]; - type = pmSourceDophotType (source); - lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); - - psLineInit (line); - psLineAdd (line, "%6.1f ", PAR[2]); - psLineAdd (line, "%6.1f ", PAR[3]); - psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->fitMag + ZERO_POINT)); - psLineAdd (line, "%03d ", PS_MIN (999, (int)(1000*dmag))); - psLineAdd (line, "%2d ", type); - psLineAdd (line, "%3.1f ", lsky); - psLineAdd (line, "%6.3f ", 99.999); // should be 'Mgal - psLineAdd (line, "%6.3f ", PS_MIN (99.999, source->apMag + ZERO_POINT)); - psLineAdd (line, "%6.2f ", PAR[4]); // should be 'FHWM x' - psLineAdd (line, "%6.2f ", PAR[5]); // should be 'FHWM y' - psLineAdd (line, "%5.1f\n", 0); // should be theta - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - return true; -} - -// elixir-style FITS table output (header + table in 1st extension) -// this format consists of a header derived from the image header -// followed by a zero-size matrix, followed by the table data -bool pmSourcesWriteCMF (psArray *sources, char *filename, psMetadata *header) { - - psArray *table; - psMetadataItem *mdi; - psMetadata *row; - psMetadata *theader; - int i, type; - psF32 *PAR, *dPAR; - float dmag, lsky; - bool status; - - // find config information for output header - float ZERO_POINT = psMetadataLookupF32 (&status, header, "ZERO_PT"); - - table = psArrayAlloc (sources->n); - table->n = 0; - - for (i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmModelSelect (source); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - dmag = dPAR[1] / PAR[1]; - type = pmSourceDophotType (source); - lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]); - - row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "X_PIX", PS_DATA_F32, "", PAR[2]); - psMetadataAdd (row, PS_LIST_TAIL, "Y_PIX", PS_DATA_F32, "", PAR[3]); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_RAW", PS_DATA_F32, "", source->fitMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_ERR", PS_DATA_F32, "", (int)(1000*dmag)); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_GAL", PS_DATA_F32, "", 32.0); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_AP", PS_DATA_F32, "", source->apMag + ZERO_POINT); - psMetadataAdd (row, PS_LIST_TAIL, "LOG_SKY", PS_DATA_F32, "", lsky); - psMetadataAdd (row, PS_LIST_TAIL, "FWHM_X", PS_DATA_F32, "", type); - psMetadataAdd (row, PS_LIST_TAIL, "FWHM_Y", PS_DATA_F32, "", PAR[4]); - psMetadataAdd (row, PS_LIST_TAIL, "THETA", PS_DATA_F32, "", PAR[5]); - psMetadataAdd (row, PS_LIST_TAIL, "DOPHOT", PS_DATA_STRING, "", "0"); - psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "123"); - - psArrayAdd (table, 100, row); - psFree (row); - } - - // write the correct number of sources in the header - psMetadataAdd (header, PS_LIST_TAIL, "NSTARS", PS_DATA_S32 | PS_META_REPLACE, "NUMBER OF STARS", table->n); - psMetadataAdd (header, PS_LIST_TAIL, "EXTEND", PS_DATA_BOOL | PS_META_REPLACE, "this file has extensions", true); - - // set NAXIS to 0 (we don't write out the data array) - mdi = psMetadataLookup (header, "NAXIS"); - mdi->data.S32 = 0; - mdi->type = PS_DATA_S32; - - // create the basic table header - theader = psMetadataAlloc (); - psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE"); - - // write out the header and table - psFits *fits = psFitsOpen (filename, "w"); - psFitsWriteHeader (header, fits); - psFitsWriteTable (fits, theader, table); - psFitsClose (fits); - - return true; -} - -/***** Text Output Methods *****/ - -bool pmSourcesWriteText (psArray *sources, char *filename) { - - char *name = (char *) psAlloc (strlen(filename) + 10); - - sprintf (name, "%s.psf.dat", filename); - pmModelWritePSFs (sources, name); - - sprintf (name, "%s.ext.dat", filename); - pmModelWriteEXTs (sources, name); - - sprintf (name, "%s.nul.dat", filename); - pmModelWriteNULLs (sources, name); - - sprintf (name, "%s.mnt.dat", filename); - pmMomentsWriteText (sources, name); - - psFree (name); - return true; -} - -// write the PSF sources to an output file -bool pmModelWritePSFs (psArray *sources, char *filename) { - - double dPos, dMag; - int i, j; - FILE *f; - psF32 *PAR, *dPAR; - pmModel *model; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("pmModelWritePSFs", 3, "can't open output file for moments%s\n", filename); - return false; - } - - // write sources with models first - for (i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - if (source->type != PM_SOURCE_STAR) continue; - model = source->modelPSF; - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - // dPos is positional error, dMag is mag error - dPos = hypot (dPAR[2], dPAR[3]); - dMag = dPAR[1] / PAR[1]; - - fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", - PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); - - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", PAR[j]); - } - fprintf (f, " : "); - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", dPAR[j]); - } - fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", - source[0].apMag, source[0].type, source[0].mode, - log10(model[0].chisq/model[0].nDOF), - source[0].moments->SN, - model[0].radius, - model[0].nDOF, - model[0].nIter); - } - fclose (f); - return true; -} - -// dump the sources to an output file -bool pmModelWriteEXTs (psArray *sources, char *filename) { - - double dPos, dMag; - int i, j; - FILE *f; - psF32 *PAR, *dPAR; - pmModel *model; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("pmModelWriteEXTs", 3, "can't open output file for moments%s\n", filename); - return false; - } - - // write sources with models first - for (i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - - if (source->type != PM_SOURCE_EXTENDED) continue; - model = pmSourceMagnitudes (source, NULL, 0.0); - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - // dPos is shape error - // XXX these are hardwired for SGAUSS - dPos = hypot ((dPAR[4] / PAR[4]), (dPAR[5] / PAR[5])); - dMag = dPAR[1] / PAR[1]; - - fprintf (f, "%7.1f %7.1f %7.1f %8.4f %7.4f %7.4f ", - PAR[2], PAR[3], PAR[0], source->fitMag, dMag, dPos); - - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", PAR[j]); - } - fprintf (f, " : "); - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", dPAR[j]); - } - fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", - source->apMag, - source[0].type, source[0].mode, - log10(model[0].chisq/model[0].nDOF), - source[0].moments->SN, - model[0].radius, - model[0].nDOF, - model[0].nIter); - } - fclose (f); - return true; -} - -// dump the sources to an output file -bool pmModelWriteNULLs (psArray *sources, char *filename) -{ - - int i; - FILE *f; - pmMoments *moment = NULL; - pmSource *source = NULL; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); - return false; - } - - pmMoments *empty = pmMomentsAlloc (); - - // write sources with models first - for (i = 0; i < sources->n; i++) { - source = sources->data[i]; - - // skip these sources (in PSF or EXT) - if (source->type == PM_SOURCE_STAR) continue; - if (source->type == PM_SOURCE_EXTENDED) continue; - - if (source->moments == NULL) { - moment = empty; - } else { - moment = source->moments; - } - - fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %8.1f %7.1f %7.1f %7.1f %4d %2d\n", - source->peak->x, source->peak->y, source->peak->counts, - source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, - source->moments->Sum, source->moments->Peak, - source->moments->Sky, source->moments->SN, - source->moments->nPixels, source->type); - } - fclose (f); - psFree (empty); - return true; -} - -// write the moments to an output file -bool pmMomentsWriteText (psArray *sources, char *filename) -{ - - int i; - FILE *f; - pmSource *source = NULL; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("pmMomentsWriteText", 3, "can't open output file for moments%s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - source = sources->data[i]; - if (source->moments == NULL) - continue; - fprintf (f, "%5d %5d %7.1f %7.1f %7.1f %6.3f %6.3f %10.1f %7.1f %7.1f %7.1f %4d %2d %#5x\n", - source->peak->x, source->peak->y, source->peak->counts, - source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, - source->moments->Sum, source->moments->Peak, - source->moments->Sky, source->moments->SN, - source->moments->nPixels, source->type, source->mode); - } - fclose (f); - return true; -} - -// write the peaks to an output file -bool pmPeaksWriteText (psArray *peaks, char *filename) { - - int i; - FILE *f; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("pmPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); - return false; - } - - for (i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - if (peak == NULL) continue; - fprintf (f, "%5d %5d %7.1f\n", - peak->x, peak->y, peak->counts); - } - fclose (f); - return true; -} - -// translations between psphot object types and dophot object types -int pmSourceDophotType (pmSource *source) { - - switch (source->type) { - - case PM_SOURCE_DEFECT: - case PM_SOURCE_SATURATED: - return (8); - - case PM_SOURCE_STAR: - if (source->mode & PM_SOURCE_SATSTAR) return (10); - if (source->mode & PM_SOURCE_POOR) return (7); - if (source->mode & PM_SOURCE_FAIL) return (4); - return (1); - - case PM_SOURCE_EXTENDED: - return (2); - - default: - return (0); - } - return (0); } Index: /trunk/psphot/src/psphotParseCamera.c =================================================================== --- /trunk/psphot/src/psphotParseCamera.c (revision 6570) +++ /trunk/psphot/src/psphotParseCamera.c (revision 6571) @@ -8,60 +8,74 @@ } } -// 2006.02.07 : no leaks! -ppFile *psphotParseCamera (ppConfig *config) { +bool psphotParseCamera (pmConfig *config) { - ppFile *input = ppFileAlloc (); + // psphot is supplied with a list of input images (may be only one image) + psArray *infiles = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); + if (infiles->n < 1) psAbort (__func__, "empty input list"); - input->filename = psMetadataLookupStr(NULL, config->arguments, "INPUT_FILE"); - psMemCopy (input->filename); // keep for external use - - // Open the input image - psLogMsg("psphot", PS_LOG_INFO, "Opening input image: %s\n", input->filename); - input->fits = psFitsOpen (input->filename, "r"); // File handle for FITS file - if (! input->fits) { - psErrorStackPrint(stderr, "Can't open input image: %s\n", input->filename); - exit(EXIT_FAILURE); + // if no camera has been specified, use the first image as a template for the rest. + if (config->camera == NULL) { + psFits *fits = psFitsOpen (infiles->data[0], "r"); + psMetadata *phu = psFitsReadHeader (NULL, fits); + config->camera = pmConfigCameraFromHeader (config->site, phu); + psFitsClose (fits); + psFree (phu); } - input->phu = psFitsReadHeader(NULL, input->fits); // FITS header - - // Get camera configuration from header if not already defined + // There's no point in continuing if we can't recognize what we've got if (! config->camera) { - config->camera = pmConfigCameraFromHeader(config->site, input->phu); - if (! config->camera) { - // There's no point in continuing if we can't recognize what we've got - psErrorStackPrint(stderr, "Can't find camera configuration!\n"); - exit(EXIT_FAILURE); - } - } else if (! pmConfigValidateCamera(config->camera, input->phu)) { - // There's no point in continuing if what we've got doesn't match what we've been told - psError(PS_ERR_IO, true, "%s does not seem to be from the specified camera.\n", - input->filename); - exit(EXIT_FAILURE); + psErrorStackPrint(stderr, "Can't find camera configuration!\n"); + exit(EXIT_FAILURE); } - // Determine the correct recipe to use (from camera or from user) - // if config->recipe is not NULL, it is a user-supplied recipe - if (!config->recipe) { - config->recipe = pmConfigRecipeFromCamera(config->camera, PSPHOT_RECIPE); - if (config->recipe == NULL) { - psErrorStackPrint(stderr, "Can't find recipe configuration!\n"); - exit(EXIT_FAILURE); - } + // XXX place the "data" element on the pmConfig structure? + config->files = psMetadataAlloc (); + + // build the template fpa, set up the basic view + pmFPA *input = pmFPAConstruct (camera); + + // assign the I/O files (potentially) needed by psphot + pmFPAfile *file = pmFPAfileDefine (config->files, camera, input, "PSPHOT.INPUT"); + + pmFPAview *view = pmFPAviewAlloc (input, config->camera, 0); + for (int i = 0; i < infiles->n; i++) { + // XXX save the phu from above? + psFits *fits = psFitsOpen (infiles->data[i], "r"); + psMetadata *phu = psFitsReadHeader (NULL, fits); + pmConfigValidateCamera (config->camera, phu); + + // set the view to the corresponding entry for this phu + // XXX should I save the phu data at the appropriate location in the fpa? + pmFPAsetView (view, phu); + + // XXX is this the correct psMD to save the filename? + name = pmFPAviewNameFromRule (file->filextra, view); + psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]); } + pmFPAfileDefine (config->files, camera, input, "PSPHOT.OUTPUT"); + pmFPAfileDefine (config->files, camera, input, "PSPHOT.RESID"); + pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_INPUT"); + pmFPAfileDefine (config->files, camera, input, "PSPHOT.PSF_OUTPUT"); + + // build the template fpa, set up the basic view + // supply the backgnd with a different camera? + pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKSUB"); + pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKGND"); + pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.BACKMDL"); + pmFPAfileDefine (config->files, camera, NULL, "PSPHOT.PSF_SAMPLE"); // recipe override values (command-line options): - psMetadataItem *item = NULL; - psMetadataIterator *iter = psMetadataIteratorAlloc (config->options, PS_LIST_HEAD, NULL); - while ((item = psMetadataGetAndIncrement (iter)) != NULL) { - psMetadataAddItem (config->recipe, item, PS_LIST_TAIL, PS_META_REPLACE); + psMetadata *options = psMetadataLookupPtr (status, config->arguments, "PSPHOT.OPTIONS"); + psMetadata *recipe = psMetadataLookupPtr (status, config->recipes, "PSPHOT"); + psMetadataIterator *iter = psMetadataIteratorAlloc (options, PS_LIST_HEAD, NULL); + while ((psMetadataItem *item = psMetadataGetAndIncrement (iter)) != NULL) { + psMetadataAddItem (recipe, item, PS_LIST_TAIL, PS_META_REPLACE); } psFree (iter); // set default recipe values here - METADATA_ADD_DEFAULT (config->recipe, Str, "FITMODE", "NONE", ""); - METADATA_ADD_DEFAULT (config->recipe, Str, "PHOTCODE", "NONE", ""); - METADATA_ADD_DEFAULT (config->recipe, Str, "BREAK_POINT", "NONE", ""); - METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_FORMAT", "CMP", ""); - METADATA_ADD_DEFAULT (config->recipe, Str, "OUTPUT_MODE", "SPLIT", ""); + // XXX this needs re-thinking... + METADATA_ADD_DEFAULT (recipe, Str, "FITMODE", "NONE", ""); + METADATA_ADD_DEFAULT (recipe, Str, "PHOTCODE", "NONE", ""); + METADATA_ADD_DEFAULT (recipe, Str, "BREAK_POINT", "NONE", ""); // Chip selection: if we are using a single chip, select it for each FPA @@ -88,9 +102,5 @@ } } - - // Construct cameras in preparation for reading - input->fpa = pmFPAConstruct(config->camera); - psTrace(__func__, 1, "Done with psphotParseCamera...\n"); - return input; + return true; } Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 6570) +++ /trunk/psphot/src/psphotReadout.c (revision 6571) @@ -2,40 +2,42 @@ // XXX 2006.02.07 : no leaks! -bool psphotReadout (pmReadout *readout, psMetadata *config) { +// XXX change 'config' to 'recipe' eventually +bool psphotReadout (pmConfig *config, pmFPAview *view) { psArray *sources = NULL; psArray *peaks = NULL; - psImage *skymodel = NULL; pmPSF *psf = NULL; bool status; + // find the currently selected readout + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); + pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + pmReadout *readout = pmFPAviewThisReadout (view, input); + // generate a background model (median, smoothed image) - skymodel = psphotImageMedian (readout, config); - - // optional output of background-subtracted image - psphot + psphotImageMedian (config, view); // find the peaks in the image - peaks = psphotFindPeaks (readout, config); + peaks = psphotFindPeaks (readout, recipe); // construct sources and measure basic stats // limit moments analysis by S/N? - sources = psphotSourceStats (readout, config, peaks); + sources = psphotSourceStats (readout, recipe, peaks); // classify sources based on moments, brightness // faint sources not classified? - psphotRoughClass (sources, config); + psphotRoughClass (sources, recipe); // mark blended peaks PS_SOURCE_BLEND - psphotBasicDeblend (sources, config); + psphotBasicDeblend (sources, recipe); // use bright stellar objects to measure PSF - psf = psphotChoosePSF (config, sources); + psf = psphotChoosePSF (sources, recipe); // linear PSF fit to peaks - psphotEnsemblePSF (readout, config, sources, psf, FALSE); + psphotEnsemblePSF (readout, sources, recipe, psf, FALSE); // non-linear PSF and EXT fit to brighter sources - psphotBlendFit (readout, config, sources, psf); + psphotBlendFit (readout, sources, recipe, psf); // replace fitted sources @@ -46,13 +48,14 @@ // linear PSF fit to remaining peaks - psphotEnsemblePSF (readout, config, sources, psf, TRUE); + psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); // measure aperture photometry corrections - psphotApResid (readout, sources, config, psf); + psphotApResid (readout, sources, recipe, psf); // calculate source magnitudes - psphotMagnitudes (config, sources, psf); + psphotMagnitudes (sources, recipe, psf); // update the header with stats results + // XXX need to do this conditionally? psMetadata *header = pmReadoutGetHeader (readout); psphotUpdateHeader (header, config); @@ -60,5 +63,5 @@ // XXX make this an option? // replace background in residual image - psphotSkyReplace (readout, skymodel); + psphotSkyReplace (config, view); // need to do something with the sources, psf, and sky @@ -67,5 +70,4 @@ // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.MEAN", PS_DATA_F32, "psphot sky mean", sky->sampleMean); // XXX : replace? status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SKY.SIGMA", PS_DATA_F32, "psphot sky stdev", sky->sampleStdev); - // XXX : what should I do with the skymodel image? // free up the local copies of the data @@ -73,5 +75,4 @@ psFree (sources); psFree (peaks); - psFree (skymodel); return true; } Index: unk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 6570) +++ (revision ) @@ -1,117 +1,0 @@ -# include "psphot.h" - -// load the image -// load or construct the mask -// apply X/Y border info -// apply the SATURATE mask? -// load or construct the weight - -pmReadout *psphotSetup (psMetadata *config, psMetadata **header) -{ - psMetadata *myHeader; - psImage *image = NULL; - psImage *weight = NULL; - psImage *mask = NULL; - psRegion region = {0,0,0,0}; // a region representing the entire array - bool status; - - psTimerStart ("psphot"); - - // load implementation-specific models - psphotModelGroupInit (); - - // setup header template, specify COMMENT and HISTORY as MULTI - // XXX EAM - shouldn't this be part of the basic FITS concept? - myHeader = psMetadataAlloc (); - psMetadataAdd (myHeader, PS_LIST_HEAD, "COMMENT", PS_DATA_METADATA_MULTI, "folder for comment", NULL); - psMetadataAdd (myHeader, PS_LIST_HEAD, "HISTORY", PS_DATA_METADATA_MULTI, "folder for history", NULL); - - // read header and image data from INPUT - char *input = psMetadataLookupPtr (&status, config, "IMAGE"); - if (!status) psAbort ("psphot", "input image not specified"); - psFits *file = psFitsOpen (input, "r"); - myHeader = psFitsReadHeader (myHeader, file); - psImage *tmpimage = psFitsReadImage (NULL, file, region, 0); - image = psImageCopy (NULL, tmpimage, PS_TYPE_F32); - psFitsClose (file); - - // psFree (input); - // XXX EAM - does the current API expect to get a copy or a view? - - // grab these values from the approrpiate location (image header if necessary) - float RDNOISE = pmConfigLookupF32 (&status, config, myHeader, "RDNOISE"); - float GAIN = pmConfigLookupF32 (&status, config, myHeader, "GAIN"); - psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_DATA_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); - psMetadataAdd (config, PS_LIST_TAIL, "GAIN", PS_DATA_F32 | PS_META_REPLACE, "gain value used", GAIN); - - // load the weight image if it is supplied, otherwise build from input - // by 'weight' we really mean variance image - char *weightName = psMetadataLookupPtr (&status, config, "WEIGHT"); - if (status == true) { - file = psFitsOpen (weightName, "r"); - weight = psFitsReadImage (NULL, file, region, 0); - psFitsClose (file); - // psFree (weightName); XXX - see psFree (input) - } else { - // build the weight image from the input image, RDNOISE, and GAIN - weight = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32); - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - weight->data.F32[iy][ix] = PS_MAX (image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN), 0.0); - } - } - } - - // load the mask if specified, otherwise construct - char *maskName = psMetadataLookupPtr (&status, config, "MASK"); - if (status == true) { - // XXX EAM require / convert mask to psU8? - file = psFitsOpen (maskName, "r"); - mask = psFitsReadImage (NULL, file, region, 0); - psFitsClose (file); - // psFree (maskName); XXX - see psFree (input) - } else { - mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); - } - - // mask the excluded outer pixels - float XMIN = psMetadataLookupF32 (&status, config, "XMIN"); - float XMAX = psMetadataLookupF32 (&status, config, "XMAX"); - float YMIN = psMetadataLookupF32 (&status, config, "YMIN"); - float YMAX = psMetadataLookupF32 (&status, config, "YMAX"); - psRegion keep = psRegionSet (XMIN, XMAX, YMIN, YMAX); - keep = psRegionForImage (image, keep); - psImageKeepRegion (mask, keep, "OR", PSPHOT_MASK_INVALID); - - // mask the saturated pixels - // XXX EAM does the mask need to grow? - float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); - for (int i = 0; status && (i < image->numRows); i++) { - for (int j = 0; j < image->numCols; j++) { - if (image->data.F32[i][j] >= SATURATION) { - mask->data.U8[i][j] |= PSPHOT_MASK_SATURATED; - } - } - } - - // mask the pixels below min threshold - float MIN_VALID = psMetadataLookupF32 (&status, config, "MIN_VALID_PIXEL"); - for (int i = 0; status && (i < image->numRows); i++) { - for (int j = 0; j < image->numCols; j++) { - if (image->data.F32[i][j] < MIN_VALID) { - mask->data.U8[i][j] |= PSPHOT_MASK_INVALID; - } - } - } - - psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); - - // return image data - pmReadout *readout = pmReadoutAlloc(NULL); - readout->image = image; - readout->weight = weight; - readout->mask = mask; - - *header = myHeader; - return (readout); -} Index: /trunk/psphot/src/psphotSkyReplace.c =================================================================== --- /trunk/psphot/src/psphotSkyReplace.c (revision 6570) +++ /trunk/psphot/src/psphotSkyReplace.c (revision 6571) @@ -1,10 +1,21 @@ # include "psphot.h" -// generate the median in NxN boxes, clipping heavily -// linear interpolation to generate full-scale model -bool psphotSkyReplace (pmReadout *readout, psImage *background) { +// in order to successfully replace the sky, we must define a corresponding file... +bool psphotSkyReplace (psMetadata *config, pmFPAview *view) { + // find the currently selected readout + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, "PSPHOT"); + pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + pmReadout *readout = pmFPAviewThisReadout (view, input); + + // select the corresponding images psImage *image = readout->image; psImage *mask = readout->mask; + + // select background pixels, from output background file, or create + background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 0, 0, PS_TYPE_F32); + if (background == NULL) { + return false; + } // replace the background model @@ -17,4 +28,5 @@ } + psFree (background); return true; } Index: /trunk/psphot/src/testPSphotLoop.c =================================================================== --- /trunk/psphot/src/testPSphotLoop.c (revision 6570) +++ /trunk/psphot/src/testPSphotLoop.c (revision 6571) @@ -1,45 +1,96 @@ # include "psphot.h" -bool psphotLoop (psphotData *data, ppConfig *config) { +bool psphotImageLoop (pmConfig *config) { - psRegion region = {0,0,0,0}; - pmFPA *fpa = data->input->fpa; - pmFPAiterator *fpi = pmFPAiteratorAlloc (fpa, region); + psMetadata *recipe = psMetadataLookupPts (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psErrorStackPrint(stderr, "Can't find recipe configuration!\n"); + exit(EXIT_FAILURE); + } + + pmFPA *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + pmFPAview *view = pmFPAviewAlloc (0); // files associated with the science image - pmFileDefine (fpi, config->camera, "PPIMAGE.INPUT"); - pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.HEAD"); - pmFileDefine (fpi, config->camera, "PPIMAGE.OUTPUT.DATA"); - pmFileDefine (fpi, config->camera, "PPIMAGE.BACKGND"); - pmFileDefine (fpi, config->camera, "PPIMAGE.BACKSUB"); - pmFileDefine (fpi, config->camera, "PPIMAGE.RESID"); - pmFileDefine (fpi, config->camera, "PPIMAGE.PSF"); - pmFileIOChecks (fpi); + pmFPAfileReadChecks (config->files, view); - while ((chip = pmChipNext (fpi)) != NULL) { + while ((chip = pmChipNext (view, input)) != NULL) { psLogMsg ("psphot", 4, "Chip %d: %x %x\n", i, chip->exists, chip->process); if (! chip->process) { continue; } - pmFileIOChecks (fpi); + pmFPAfileReadChecks (config->files, view); - while ((cell = pmCellNext (fpi)) != NULL) { + while ((cell = pmCellNext (view, input)) != NULL) { psLogMsg ("psphot", 4, "Cell %d: %x %x\n", j, cell->exists, cell->process); if (! cell->process) { continue; } - pmFileIOChecks (fpi); + pmFPAfileReadChecks (config->files, view); + + // XXX optional mask and weight input image should be loaded here? + // this sets the weight map and basic mask applying CELL.BAD and CELL.SATURATION + pmCellSetWeights(cell); + + // I have a valid mask, now mask in the analysis region of interest + pmCellSetMask (cell, recipe); // process each of the readouts - while ((readout = pmReadoutNext (fpi)) != NULL) { - pmFileIOChecks (fpi); + while ((readout = pmReadoutNext (view, input)) != NULL) { + pmFPAfileReadChecks (config->files, view); + // run a single-model test if desired + // XXX move this to psphotReadout?? + // psphotModelTest (readout, recipe); + // run the actual photometry analysis - psphotReadout (readout, config); + psphotReadout (config, view); + + pmFPAfileWriteChecks (config->files, view); // write out the desired output dataset(s) - psphotOutput (readout, config); + // psphotOutput (view, recipe); } + pmFPAfileWriteChecks (config->files, view); } + pmFPAfileWriteChecks (config->files, view); } - pmFileClose (fpi); + pmFPAfileWriteChecks (config->files, view); return true; } + +// I/O files related to psphot: +// PSPHOT.INPUT : input image file(s) +// PSPHOT.RESID : residual image +// PSPHOT.OUTPUT : output object tables (object) + +// PSPHOT.BACKSUB : background subtracted image +// PSPHOT.BACKGND : background model (full-scale image?) +// PSPHOT.BACKMDL : background model (binned image?) +// PSPHOT.PSF : sample PSF images + + +/** + +filename | pmFPAfile | pmFPA | pmFPAview +chip00.f | PSPHOT.IN | input | view +chip01.f | PSPHOT.IN | input | view +chip02.f | PSPHOT.IN | input | view +chip03.f | PSPHOT.IN | input | view + +out00.f | PSPHOT.OUT | input | view +out01.f | PSPHOT.OUT | input | view +out02.f | PSPHOT.OUT | input | view +out03.f | PSPHOT.OUT | input | view + +obj00.f | PSPHOT.OBJ | input | view +obj01.f | PSPHOT.OBJ | input | view +obj02.f | PSPHOT.OBJ | input | view +obj03.f | PSPHOT.OBJ | input | view + +bin00.f | PSPHOT.BIN | binned | view +bin01.f | PSPHOT.BIN | binned | view +bin02.f | PSPHOT.BIN | binned | view +bin03.f | PSPHOT.BIN | binned | view + +mosaic.f | PSPHOT.MOS | mosaic | view + +**/