Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 4945) +++ /trunk/psphot/Makefile (revision 4946) @@ -3,5 +3,6 @@ @echo "USAGE: make psphot" -CC = gcc -g -Wall -O2 +# CC = gcc -g -Wall -O2 +CC = gcc -g -Wall SRC = src BIN = bin @@ -18,25 +19,25 @@ PSPHOT = \ -$(SRC)/psphot.$(ARCH).o \ -$(SRC)/pspsf.$(ARCH).o \ -$(SRC)/psphot-utils.$(ARCH).o \ -$(SRC)/psPolynomials.$(ARCH).o \ -$(SRC)/psUtils.$(ARCH).o \ -$(SRC)/load_args.$(ARCH).o \ -$(SRC)/setup.$(ARCH).o \ -$(SRC)/image_stats.$(ARCH).o \ -$(SRC)/find_peaks.$(ARCH).o \ -$(SRC)/source_moments.$(ARCH).o \ -$(SRC)/choose_psf_model.$(ARCH).o \ -$(SRC)/apply_psf_model.$(ARCH).o \ -$(SRC)/test_psf_scatter.$(ARCH).o \ -$(SRC)/mark_psf_source.$(ARCH).o \ -$(SRC)/subtract_psf_source.$(ARCH).o \ -$(SRC)/by_SN.$(ARCH).o \ -$(SRC)/fit_galaxies.$(ARCH).o \ -$(SRC)/LocalSky.$(ARCH).o \ -$(SRC)/basic_classes.$(ARCH).o \ -$(SRC)/psLine.$(ARCH).o \ -$(SRC)/output.$(ARCH).o +$(SRC)/psphot.$(ARCH).o \ +$(SRC)/psphotArguments.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psphotImageStats.$(ARCH).o \ +$(SRC)/psphotSourceStats.$(ARCH).o \ +$(SRC)/psphotChoosePSF.$(ARCH).o \ +$(SRC)/psphotApplyPSF.$(ARCH).o \ +$(SRC)/psphotFitGalaxies.$(ARCH).o \ +$(SRC)/psphotOutput.$(ARCH).o \ +$(SRC)/psphotMarkPSF.$(ARCH).o \ +$(SRC)/psphotSubtractPSF.$(ARCH).o \ +$(SRC)/pspsf.$(ARCH).o \ +$(SRC)/psLibUtils.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ +$(SRC)/pmSourceUtils.$(ARCH).o \ +$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \ +$(SRC)/psLine.$(ARCH).o \ +$(SRC)/psPolynomials.$(ARCH).o \ +$(SRC)/psMinimize.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psphotSortBySN.$(ARCH).o FITSOURCE = \ Index: unk/psphot/src/LocalSky.c =================================================================== --- /trunk/psphot/src/LocalSky.c (revision 4945) +++ (revision ) @@ -1,590 +1,0 @@ -# include "psphot.h" - -bool pmSourceDefinePixels(psSource *mySource, - const psImageData *imdata, - psF32 x, - psF32 y, - psF32 Radius) -{ - psRegion *srcRegion; - - // Grab a subimage of the original image of size (2 * outerRadius). - srcRegion = psRegionSquare (x, y, Radius); - srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion); - - // use these when psImageSubset is updated - // psImage *subImage = psImageSubset(imdata->image, srcRegion); - // psImage *subImageMask = psImageSubset(imdata->mask, srcRegion); - // psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion); - - mySource->pixels = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - mySource->noise = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - mySource->mask = psImageSubset(imdata->mask, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); - - psFree (srcRegion); - return(mySource); -} - -// this sets and clears bit 0x80 -bool pmSourceLocalSky_EAM (psSource *source, - psStatsOptions statsOptions, - psF32 Radius) -{ - - psImage *image = source->pixels; - psImage *mask = source->mask; - psPeak *peak = source->peak; - psRegion *srcRegion; - - srcRegion = psRegionSquare (peak->x, peak->y, Radius); - srcRegion = psRegionForImage (srcRegion, mask, srcRegion); - psImageMaskRegion (mask, srcRegion, OR, 0x80); - - psStats *myStats = psStatsAlloc(statsOptions); - myStats = psImageStats(myStats, image, mask, 0xff); - psImageMaskRegion (mask, srcRegion, AND, 0x7f); - - psF64 tmpF64; - p_psGetStatValue(myStats, &tmpF64); - - psFree (myStats); - psFree (srcRegion); - - if (isnan(tmpF64)) return (false); - source->moments = pmMomentsAlloc(); - source->moments->Sky = (psF32) tmpF64; - return (true); -} - -# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 -# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 -bool pmSourceFitModel_EAM(psSource *source, - psModel *model, - const bool PSF) -{ - PS_PTR_CHECK_NULL(source, false); - PS_PTR_CHECK_NULL(source->moments, false); - PS_PTR_CHECK_NULL(source->peak, false); - PS_PTR_CHECK_NULL(source->pixels, false); - PS_PTR_CHECK_NULL(source->mask, false); - PS_PTR_CHECK_NULL(source->noise, false); - psBool fitStatus = true; - psBool onPic = true; - psBool rc = true; - psF32 Ro, ymodel; - - - // XXX EAM : is it necessary for the mask & noise to exist? the - // tests below could be conditions (!NULL) - - psModelFunc modelFunc = psModelFunc_GetFunction (model->type); - psModelLimits modelLimits = psModelLimits_GetFunction (model->type); - - psVector *params = model->params; - psVector *dparams = model->dparams; - psVector *paramMask = NULL; - - psVector *beta_lim = NULL; - psVector *params_min = NULL; - psVector *params_max = NULL; - - int nParams = PSF ? params->n - 4 : params->n; - psF32 So = params->data.F32[0]; - - // find the number of valid pixels - // XXX EAM : this loop and the loop below could just be one pass - // using the psArrayAdd and psVectorExtend functions - psS32 count = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - count++; - } - } - } - if (count < nParams + 1) { - psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); - return(false); - } - - // construct the coordinate and value entries - psArray *x = psArrayAlloc(count); - psVector *y = psVectorAlloc(count, PS_TYPE_F32); - psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); - psS32 tmpCnt = 0; - for (psS32 i = 0; i < source->pixels->numRows; i++) { - for (psS32 j = 0; j < source->pixels->numCols; j++) { - if (source->mask->data.U8[i][j] == 0) { - psVector *coord = psVectorAlloc(2, PS_TYPE_F32); - // XXX: Convert i/j to image space: - // XXX EAM: coord order is (x,y) == (col,row) - coord->data.F32[0] = (psF32) (j + source->pixels->col0); - coord->data.F32[1] = (psF32) (i + source->pixels->row0); - x->data[tmpCnt] = (psPtr *) coord; - y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; - - ymodel = modelFunc (NULL, model->params, coord); - - // this test enhances the noise based on deviation from the model flux - Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); - yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro); - // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then - // the minimization function calculates sq() - tmpCnt++; - } - } - } - - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - - // PSF model only fits first 4 parameters, FLT model fits all - if (PSF) { - paramMask = psVectorAlloc (params->n, PS_TYPE_U8); - for (int i = 0; i < 4; i++) { - paramMask->data.U8[i] = 0; - } - for (int i = 4; i < paramMask->n; i++) { - paramMask->data.U8[i] = 1; - } - } - - psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); - modelLimits (&beta_lim, ¶ms_min, ¶ms_max); - for (int i = 0; i < params->n; i++) { - covar->data.F64[0][i] = beta_lim->data.F32[i]; - covar->data.F64[1][i] = params_min->data.F32[i]; - covar->data.F64[2][i] = params_max->data.F32[i]; - } - - psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); - fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if ((paramMask != NULL) && paramMask->data.U8[i]) continue; - dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); - } - - // XXX EAM: we need to do something (give an error?) if rc is false - // XXX EAM: psMinimizeLMChi2 does not check convergence - - // XXX models can go insane: reject these - onPic &= (params->data.F32[2] >= source->pixels->col0); - onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); - onPic &= (params->data.F32[3] >= source->pixels->row0); - onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); - - // XXX EAM: save the resulting chisq, nDOF, nIter - model->chisq = myMin->value; - model->nIter = myMin->iter; - model->nDOF = y->n - nParams; - - // XXX EAM get the Gauss-Newton distance for fixed model parameters - if (paramMask != NULL) { - psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); - psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); - for (int i = 0; i < dparams->n; i++) { - if (!paramMask->data.U8[i]) continue; - dparams->data.F32[i] = delta->data.F64[i]; - } - } - - psFree(paramMask); - psFree(x); - psFree(y); - psFree(myMin); - - rc = (onPic && fitStatus); - return(rc); -} - -/****************************************************************************** -pmSourceMoments(source, radius): this function takes a subImage defined in the -psSource data structure, along with the peak location, and determines the -various moments associated with that peak. - -Requires the following to have been created: - psSource - psSource->peak - psSource->pixels - -XXX: The peak calculations are done in image coords, not subImage coords. - -XXX: mask values? -*****************************************************************************/ -# define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) - -bool pmSourceMoments_EAM(psSource *source, - psF32 radius) -{ - PS_PTR_CHECK_NULL(source, NULL); - PS_PTR_CHECK_NULL(source->peak, NULL); - PS_PTR_CHECK_NULL(source->pixels, NULL); - PS_PTR_CHECK_NULL(source->noise, NULL); - PS_FLOAT_COMPARE(0.0, radius, NULL); - - // - // XXX: Verify the setting for sky if source->moments == NULL. - // - psF32 sky = 0.0; - if (source->moments == NULL) { - source->moments = pmMomentsAlloc(); - } else { - sky = source->moments->Sky; - } - - // - // Sum = SUM (z - sky) - // X1 = SUM (x - xc)*(z - sky) - // X2 = SUM (x - xc)^2 * (z - sky) - // XY = SUM (x - xc)*(y - yc)*(z - sky) - // - psF32 peakPixel = -PS_MAX_F32; - psS32 numPixels = 0; - psF32 Sum = 0.0; - psF32 X1 = 0.0; - psF32 Y1 = 0.0; - psF32 X2 = 0.0; - psF32 Y2 = 0.0; - psF32 XY = 0.0; - psF32 x = 0; - psF32 y = 0; - psF32 R2 = PS_SQR(radius); - - psF32 xPeak = source->peak->x; - psF32 yPeak = source->peak->y; - - // XXX why do I get different results for these two methods of finding Sx? - // XXX Sx, Sy would be better measured if we clip pixels close to sky - // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? - // We loop through all pixels in this subimage (source->pixels), and for each - // pixel that is not masked, AND within the radius of the peak pixel, we - // proceed with the moments calculation. need to do two loops for a - // numerically stable result. first loop: get the sums. - // XXX EAM : mask == 0 is valid - - for (psS32 row = 0; row < source->pixels->numRows ; row++) { - for (psS32 col = 0; col < source->pixels->numCols ; col++) { - if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue; - - psF32 xDiff = col + source->pixels->col0 - xPeak; - psF32 yDiff = row + source->pixels->row0 - yPeak; - - if (!VALID_RADIUS(xDiff, yDiff)) continue; - - psF32 pDiff = source->pixels->data.F32[row][col] - sky; - - // XXX EAM : check for valid S/N in pixel - if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue; - - Sum += pDiff; - X1 += xDiff * pDiff; - Y1 += yDiff * pDiff; - XY += xDiff * yDiff * pDiff; - - X2 += PS_SQR(xDiff) * pDiff; - Y2 += PS_SQR(yDiff) * pDiff; - - peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); - numPixels++; - } - } - // XXX EAM - the limit is a bit arbitrary. make it user defined? - if ((numPixels < 3) || (Sum <= 0)) { - psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); - return (false); - } - - psTrace (".psModules.pmSourceMoments", 5, - "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", - sky, Sum, X1, Y1, X2, Y2, XY, numPixels); - - // - // first moment X = X1/Sum + xc - // second moment X = sqrt (X2/Sum - (X1/Sum)^2) - // Sxy = XY / Sum - // - x = X1/Sum; - y = Y1/Sum; - if ((fabs(x) > radius) || (fabs(y) > radius)) { - psTrace (".psModules.pmSourceMoments", 5, - "large centroid swing; invalid peak %d, %d\n", - source->peak->x, source->peak->y); - return (false); - } - - source->moments->x = x + xPeak; - source->moments->y = y + yPeak; - - source->moments->Sxy = XY/Sum - x*y; - source->moments->Sum = Sum; - source->moments->Peak = peakPixel; - source->moments->nPixels = numPixels; - - // XXX EAM : these values can be negative, so we need to limit the range - source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); - source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); - - psTrace (".psModules.pmSourceMoments", 4, - "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", - sky, Sum, source->moments->x, source->moments->y, - source->moments->Sx, source->moments->Sy, source->moments->Sxy); - - return(true); -} - -bool pmModelFitStatus (psModel *model) { - - bool status; - - psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type); - status = statusFunc (model); - - return (status); -} - -// XXX EAM this implementation of MinLM includes limits on params & dparams -psBool psMinimizeLMChi2_EAM(psMinimization *min, - psImage *covar, - psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func) -{ - PS_PTR_CHECK_NULL(min, NULL); - PS_VECTOR_CHECK_NULL(params, NULL); - PS_VECTOR_CHECK_EMPTY(params, NULL); - PS_PTR_CHECK_NULL(x, NULL); - PS_VECTOR_CHECK_NULL(y, NULL); - PS_VECTOR_CHECK_EMPTY(y, NULL); - PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); - PS_PTR_CHECK_NULL(func, NULL); - - // this function has test and current values for several things - // the current best value is in lower case - // the next guess value is in upper case - - // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32); - psVector *dy = NULL; - psF64 Chisq = 0.0; - psF64 lambda = 0.001; - - psVector *beta_lim = NULL; - psVector *param_min = NULL; - psVector *param_max = NULL; - - // if we are provided a covar image, we expect to find these three vectors in first three rows - if (covar != NULL) { - beta_lim = psVectorAlloc (params->n, PS_TYPE_F32); - param_min = psVectorAlloc (params->n, PS_TYPE_F32); - param_max = psVectorAlloc (params->n, PS_TYPE_F32); - for (int i = 0; i < params->n; i++) { - beta_lim->data.F32[i] = covar->data.F64[0][i]; - param_min->data.F32[i] = covar->data.F64[1][i]; - param_max->data.F32[i] = covar->data.F64[2][i]; - } - psImageRecycle (covar, params->n, params->n, PS_TYPE_F64); - } - - - // why is this needed here??? the initial guess on params is provided by the user - Params = psVectorCopy (Params, params, PS_TYPE_F32); - - // the user provides the error or NULL. we need to convert - // to appropriate weights - dy = psVectorAlloc (y->n, PS_TYPE_F32); - if (yErr != NULL) { - for (int i = 0; i < dy->n; i++) { - if (yErr->data.F32[i] == 0.0) { - dy->data.F32[i] = 1.0; - // mask this? bad pixel, obviously... - } else { - dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]); - } - } - } else { - for (int i = 0; i < dy->n; i++) { - dy->data.F32[i] = 1.0; - } - } - - // calculate initial alpha and beta, set chisq (min->value) - min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); - if (isnan(min->value)) { - min->iter = min->maxIter; - return (false); - } - # ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), params, "params guess"); - } - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { - p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); - } - # endif /* PS_NO_TRACE */ - - // iterate until the tolerance is reached, or give up - while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - - // set a new guess for Alpha, Beta, Params - p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda); - - // measure linear model prediction - psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda); - - # ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); - } - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { - p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); - } - # endif /* PS_NO_TRACE */ - - // calculate Chisq for new guess, update Alpha & Beta - Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func); - - // XXX EAM alternate convergence criterion: - // compare the delta (min->value - Chisq) with the - // expected delta from the linear model (dLinear) - // accept new guess (if improvement), or increase lambda - psF64 rho = (min->value - Chisq) / dLinear; - - psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho); - # ifndef PS_NO_TRACE - // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { - p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); - p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); - p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); - } - # endif /* PS_NO_TRACE */ - - /* if (Chisq < min->value) { */ - if (rho > 0.0) { - min->lastDelta = (min->value - Chisq) / (dy->n - params->n); - min->value = Chisq; - alpha = psImageCopy (alpha, Alpha, PS_TYPE_F64); - beta = psVectorCopy (beta, Beta, PS_TYPE_F64); - params = psVectorCopy (params, Params, PS_TYPE_F32); - lambda *= 0.1; - } else { - lambda *= 10.0; - } - min->iter ++; - } - psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); - - // construct & return the covariance matrix (if requested) - if (covar != NULL) { - p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0); - } - - // free the internal temporary data - psFree (alpha); - psFree (Alpha); - psFree (beta); - psFree (Beta); - psFree (Params); - psFree (dy); - - if (min->iter == min->maxIter) { - return (false); - } - return (true); -} - -// XXX EAM : can we use static copies of LUv, LUm, A? -psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, - psVector *Beta, - psVector *Params, - const psImage *alpha, - const psVector *beta, - const psVector *params, - const psVector *paramMask, - const psVector *beta_lim, - const psVector *params_min, - const psVector *params_max, - psF64 lambda) -{ - - # define USE_LU_DECOMP 1 - # if (USE_LU_DECOMP) - psVector *LUv = NULL; - psImage *LUm = NULL; - psImage *A = NULL; - psF32 det; - - // LU decomposition version - psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n"); - - // set new guess values (creates matrix A) - A = psImageCopy (NULL, alpha, PS_TYPE_F64); - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); - } - - // solve A*beta = Beta (Alpha = 1/A) - // these operations do not modify the input values (creates LUm, LUv) - LUm = psMatrixLUD (NULL, &LUv, A); - Beta = psMatrixLUSolve (Beta, LUm, beta, LUv); - Alpha = psMatrixInvert (Alpha, A, &det); - - # else - // gauss-jordan version - psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version"); - - // set new guess values (creates matrix A) - Beta = psVectorCopy (Beta, beta, PS_TYPE_F64); - Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); - } - - psGaussJordan (Alpha, Beta); - # endif - - // apply Beta to get new Params values - for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; - // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; - // continue; - // compare Beta to beta limits - if (beta_lim != NULL) { - if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) { - Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]); - } - } - Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; - // compare new params to param limits - if (params_max != NULL) { - Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]); - } - if (params_min != NULL) { - Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]); - } - } - - # if (USE_LU_DECOMP) - psFree (A); - psFree (LUm); - psFree (LUv); - # endif - - return true; -} - Index: unk/psphot/src/apply_psf_model.c =================================================================== --- /trunk/psphot/src/apply_psf_model.c (revision 4945) +++ (revision ) @@ -1,93 +1,0 @@ -# include "psphot.h" - -// fit psf model to all objects -// PSFSTAR objects will be refitted -// run this function to a specific flux limit? - -bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) -{ - bool status; - float x; - float y; - int Nfit = 0; - int Nsub = 0; - int Niter = 0; - - psMemoryId id = psMemGetId (); - psTimerStart ("psphot"); - - // we may set this differently here from the value used to mark likely saturated stars - float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE"); - float FIT_NSIGMA = psMetadataLookupF32 (&status, config, "FIT_NSIGMA"); - float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING"); - float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); - float OUTER = psMetadataLookupF32 (&status, config, "OUTER_RADIUS"); - float FIT_MIN_SN = psMetadataLookupF32 (&status, config, "FIT_MIN_SN"); - float FIT_MAX_CHI = psMetadataLookupF32 (&status, config, "FIT_MAX_CHI"); - - // set the object surface-brightness limit for fitted pixels - float FLUX_LIMIT = FIT_NSIGMA * sky->sampleStdev; - psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); - - // this function specifies the radius at this the model hits the given flux - psModelRadius modelRadius = psModelRadius_GetFunction (psf->type); - - for (int i = 0; i < sources->n; i++) { - - psSource *source = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - if (source->type == PS_SOURCE_DEFECT) continue; // XX should I try these anyway? - if (source->type == PS_SOURCE_SATURATED) continue; - - // use the source moments, etc to guess basic model parameters - psModel *modelFLT = pmSourceModelGuess (source, psf->type); - - // set PSF parameters for this model - psModel *model = psModelFromPSF (modelFLT, psf); - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - psFree (modelFLT); - - // set the fit radius based on the object flux limit and the model - // FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING; - if (isnan(model->radius)) { - psAbort ("apply_psf_model", "error in radius"); - } - - // check if we need to redefine the pixels - if (model->radius > OUTER) { - // allocate image, noise, mask arrays for each peak (square of radius OUTER) - pmSourceDefinePixels (source, imdata, x, y, OUTER); - } - - // set the valid/invalid pixel mask - psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80); - status = pmSourceFitModel (source, model, true); - psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f); - if (!status || (model->params->data.F32[1] < 0)) { - psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); - source->type = PS_SOURCE_FAIL_FIT_PSF; // better choice? - psFree (model); - continue; - } - - source->modelPSF = model; - Niter += model[0].nIter; - Nfit ++; - - mark_psf_source (source, shapeNsigma, FIT_MIN_SN, FIT_MAX_CHI, SATURATE); - if (subtract_psf_source (source)) { - Nsub ++; - } - } - - // XXX this test compares dSx,dSy distribution with model. - // XXX it has some problems - // test_psf_scatter (sources); - - psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); - psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub); - return (true); -} Index: unk/psphot/src/basic_classes.c =================================================================== --- /trunk/psphot/src/basic_classes.c (revision 4945) +++ (revision ) @@ -1,13 +1,0 @@ -# include "psphot.h" - -bool basic_classes (psArray *sources, psMetadata *config) { - - pmPSFClump psfClump; - - // group into STAR, COSMIC, GALAXY, SATURATED - // XXX - this has 1 leak - psfClump = pmSourcePSFClump (sources, config); - pmSourceRoughClass (sources, config, psfClump); - - return (true); -} Index: unk/psphot/src/by_SN.c =================================================================== --- /trunk/psphot/src/by_SN.c (revision 4945) +++ (revision ) @@ -1,19 +1,0 @@ -# include "psphot.h" - -// sort by SN (descending) -int by_SN (const void **a, const void **b) -{ - psSource *A = *(psSource **)a; - psSource *B = *(psSource **)b; - - psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN; - psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN; - if (isnan (fA)) fA = 0; - if (isnan (fB)) fB = 0; - - psF32 diff = fA - fB; - if (diff > FLT_EPSILON) return (-1); - if (diff < FLT_EPSILON) return (+1); - return (0); -} - Index: unk/psphot/src/choose_psf_model.c =================================================================== --- /trunk/psphot/src/choose_psf_model.c (revision 4945) +++ (revision ) @@ -1,84 +1,0 @@ -# include "psphot.h" - -// test PSF models and select best option - -pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *skystats) -{ - bool status; - char *modelName; - char key[64]; - pmPSF_Test *test = NULL; - psArray *stars = NULL; - - // select the possible PSF stars (pointers to original sources) - stars = psArrayAlloc (sources->n); - stars->n = 0; - - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - if (source->type != PS_SOURCE_PSFSTAR) continue; - psArrayAdd (stars, 200, source); - } - psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n); - - // define model fit pixels - // check that FIT_RADIUS < OUTER_RADIUS - float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS"); - - // the user specifies 1 or more PSF model options - // run test on each model option - int Ntest = psMetadataLookupS32 (&status, config, "PSF_MODEL_N"); - psArray *tests = psArrayAlloc (Ntest); - for (int i = 0; i < tests->n; i++) { - snprintf (key, 64, "PSF_MODEL_%d", i); - modelName = psMetadataLookupPtr (&status, config, key); - - tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS); - psFree (modelName); - - // make this an option: - // sprintf (filename, "metric.%d.dat", i); - // DumpPSFTestData (tests->data[i], filename); - } - psFree (stars); - - // select the best of the test models - // here we are using the clippedStdev on the metric as the indicator - test = tests->data[0]; - int bestN = 0; - float bestM = test->dApResid; - for (int i = 1; i < Ntest; i++) { - test = tests->data[i]; - float M = test->dApResid; - if (M < bestM) { - bestM = M; - bestN = i; - } - } - // keep only the selected test: - test = tests->data[bestN]; - modelName = psModelGetType (test->modelType); - psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid); - - // XXX I am having trouble tracking down the double frees vs leaks - # if (0) - // we are keeping the modelPSF fits from the PSF test, - // but we replace these later on (apply_psf_model) - // set source->model based on best psf model fit - for (int i = 0; i < test->sources->n; i++) { - psSource *source = test->sources->data[i]; - // drop masked sources from PSFSTAR list - if (test->mask->data.U8[i]) { - source->type = PS_SOURCE_OTHER; - source->modelPSF = NULL; - } else { - // source->modelPSF = psMemCopy(test->modelPSF->data[i]); - source->modelPSF = test->modelPSF->data[i]; - } - } - # endif - - pmPSF *psf = psMemCopy(test->psf); - psFree (tests); - return (psf); -} Index: unk/psphot/src/find_peaks.c =================================================================== --- /trunk/psphot/src/find_peaks.c (revision 4945) +++ (revision ) @@ -1,43 +1,0 @@ -# include "psphot.h" - -psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky) -{ - bool status = false; - float NSIGMA; - float SIGMA; - float threshold; - - // smooth the image - - psTimerStart ("psphot"); - - SIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_SIGMA"); - NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA"); - - psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); - psImageSmooth (smooth, SIGMA, NSIGMA); - psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); - - psTimerStart ("psphot"); - - // set peak threshold - NSIGMA = psMetadataLookupF32 (&status, config, "PEAK_NSIGMA"); - threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; - psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold); - - // find the peaks in the smoothed image - psArray *peaks = pmFindImagePeaks (smooth, threshold); - if (peaks == NULL) psAbort ("find peaks", "no peaks found"); - psFree (smooth); - - // make this optional - // DumpPeaks (peaks, "peaks.dat"); - - psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); - - return (peaks); -} - -// In this function, we smooth the image, then search for the peaks -// Should we also subtract a super-binned image? (as an option?) -// We need to gracefully handle no source detections Index: unk/psphot/src/fit_galaxies.c =================================================================== --- /trunk/psphot/src/fit_galaxies.c (revision 4945) +++ (revision ) @@ -1,122 +1,0 @@ -# include "psphot.h" - -bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) -{ - bool status, goodfit; - float x; - float y; - float sky; - int Nfit = 0; - int Nfail = 0; - int Niter = 0; - - float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); - float snFaint = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM"); - float OUTER = psMetadataLookupF32 (&status, config, "OUTER_RADIUS"); - float FIT_NSIGMA = psMetadataLookupF32 (&status, config, "FIT_NSIGMA"); - float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING"); - char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); - - float FLUX_LIMIT = FIT_NSIGMA * skyStats->sampleStdev; - - psModelType modelType = psModelSetType (modelName); - psModelRadius modelRadius = psModelRadius_GetFunction (modelType); - - psTimerStart ("psphot"); - - float dt = 0; - float dT = 0; - float dT2 = 0; - int dTn = 0; - - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - - // sources which should not be fitted - // skip all valid stars - if (source->type == PS_SOURCE_PSFSTAR) continue; - if (source->type == PS_SOURCE_SATSTAR) continue; - if (source->type == PS_SOURCE_GOODSTAR) continue; - - // skip all likely defects - if (source->type == PS_SOURCE_DEFECT) continue; - if (source->type == PS_SOURCE_SATURATED) continue; - - // skip poorly fitted stars - if (source->type == PS_SOURCE_FAINTSTAR) continue; - if (source->type == PS_SOURCE_POOR_FIT_PSF) continue; - - // XXX when do we pick these up again? - if (source->moments->SN < snFaint) { - source->type = PS_SOURCE_FAINT_GALAXY; // better choice? - continue; - } - - // recalculate the source moments using the galaxy radius (larger) - status = pmSourceMoments_EAM (source, MOMENT_R); - if (!status) { - source->type = PS_SOURCE_DROP_GALAXY; // better choice? - continue; - } - - // use the source moments, etc to guess basic model parameters - psModel *model = pmSourceModelGuess (source, modelType); - - x = model->params->data.F32[2]; - y = model->params->data.F32[3]; - - // set the fit radius based on the object flux limit and the model - // FLUX_LIMIT should be set based on local sky model (not global median) - model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING; - if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); - - if (model->radius > OUTER) { - // allocate image, noise, mask arrays for each peak (square of radius OUTER) - pmSourceDefinePixels (source, imdata, x, y, model->radius); - } - - // fit as FLT, not PSF (skip poor fits) - psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80); - psTimerStart ("galfit"); - status = pmSourceFitModel_EAM (source, model, false); - dt = psTimerMark ("galfit"); - dT += dt; - dT2 += dt*dt; - dTn ++; - psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f); - if (!status) { - // if the fit fails, we need to change the classification - psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PS_SOURCE_FAIL_FIT_GAL; // better choice? - source->modelFLT = model; - Nfail ++; - continue; - } - - goodfit = pmModelFitStatus (model); - if (!goodfit) { - // if the fit fails, we need to change the classification - psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); - source->type = PS_SOURCE_POOR_FIT_GAL; // better choice? - source->modelFLT = model; - Nfail ++; - continue; - } - - source->type = PS_SOURCE_GALAXY; - source->modelFLT = model; - Niter += model[0].nIter; - Nfit++; - - // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - - } - dt = dT / dTn; - fprintf (stderr, "dt: %f +/- %f\n", dt, sqrt(dT2/dTn - PS_SQR(dt))); - psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter); - return (true); -} Index: unk/psphot/src/image_stats.c =================================================================== --- /trunk/psphot/src/image_stats.c (revision 4945) +++ (revision ) @@ -1,56 +1,0 @@ -# include "psphot.h" - -psStats *image_stats (psImageData *imdata, psMetadata *config) -{ - psStats *stats = NULL; - psStats *sky = NULL; - - psImage *image = imdata->image; - - // get image stats on a subset of the image (random 1e5 pts) - psTimerStart ("psphot"); - - // this section fixes a bug in the psImageStats definition: - // it should limit the number of elements used for the initial - // pass on the stats - { - bool status = false; - int Npixels = image->numRows*image->numCols; - int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "NSUBSET")); - psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32); - psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); - - // choose Nsubset points between 0 and Nx*Ny, convert to coords - Npixels = image->numRows*image->numCols; - for (int i = 0; i < Nsubset; i++) { - double frnd = psRandomUniform (rnd); - int pixel = Npixels * frnd; - int ix = pixel / image->numCols; - int iy = pixel % image->numCols; - subset->data.F32[i] = image->data.F32[iy][ix]; - } - - // this should use ROBUST not SAMPLE median - // robust median is broken in pslib (0.5) - stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); - stats = psVectorStats (stats, subset, NULL, NULL, 0); // use a mask? - psFree (subset); - psFree (rnd); - } - - // we store these values in mean,stdev - bool status = false; - float NOISE = psMetadataLookupF32 (&status, config, "RDNOISE"); - float GAIN = psMetadataLookupF32 (&status, config, "GAIN"); - - // convert instrumental background to poisson stats - sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - sky->sampleMean = stats->sampleMedian; - sky->sampleStdev = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE)); - - psLogMsg ("psphot", 3, "stats: %f sec\n", psTimerMark ("psphot")); - psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev); - - psFree (stats); - return (sky); -} Index: unk/psphot/src/load_args.c =================================================================== --- /trunk/psphot/src/load_args.c (revision 4945) +++ (revision ) @@ -1,64 +1,0 @@ -# include "psphot.h" - -psMetadata *load_args (int *argc, char **argv) { - - int N, Nfail; - - // basic pslib options - psLogSetFormat ("M"); - psLogArguments (argc, argv); - psTraceArguments (argc, argv); - - // identify options in args list - these go on config - - char *mask = NULL; - if ((N = get_argument (*argc, argv, "-mask"))) { - remove_argument (N, argc, argv); - mask = psStringCopy (argv[N]); - remove_argument (N, argc, argv); - } - - char *noise = NULL; - if ((N = get_argument (*argc, argv, "-noise"))) { - remove_argument (N, argc, argv); - noise = psStringCopy (argv[N]); - remove_argument (N, argc, argv); - } - - char *diff = NULL; - if ((N = get_argument (*argc, argv, "-diff"))) { - remove_argument (N, argc, argv); - diff = psStringCopy (argv[N]); - remove_argument (N, argc, argv); - } - - if (*argc != 4) usage (); - - // load config information - // XXX - 3 leaks from this - psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[3], FALSE); - - // identify input image & optional noise & mask images - psMetadataAddStr (config, PS_LIST_HEAD, "INPUT", "", argv[1]); - psMetadataAddStr (config, PS_LIST_HEAD, "OUTPUT", "", argv[2]); - if (mask != NULL) { - psMetadataAddStr (config, PS_LIST_HEAD, "MASK", "", mask); - } - if (noise != NULL) { - psMetadataAddStr (config, PS_LIST_HEAD, "NOISE", "", noise); - } - if (diff != NULL) { - psMetadataAddStr (config, PS_LIST_HEAD, "DIFF_IMAGE", "", diff); - } - return (config); -} - -int usage () { - - fprintf (stderr, "USAGE: psphot (image) (output) (config)\n"); - fprintf (stderr, "options: \n"); - fprintf (stderr, " -mask (filename)\n"); - fprintf (stderr, " -noise (filename)\n"); - fprintf (stderr, " -diff (filename)\n"); - exit (2); -} Index: unk/psphot/src/mark_psf_source.c =================================================================== --- /trunk/psphot/src/mark_psf_source.c (revision 4945) +++ (revision ) @@ -1,89 +1,0 @@ -# include "psphot.h" - -// identify objects consistent with PSF shape/magnitude distribution -// we expect dparams[4],dparams[5] to have a scatter of: -// sigma_x / (S/N) * sqrt(2) -// 1 / (params[4],params[5])*(S/N) - -// any objects which is consistent with the PSF should have -// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) -// this should include a minimum dsx buffer for the brighter objects -// saturated stars should fall outside (but are already IDed) -// galaxies should be larger, cosmic rays smaller, but need to test? - -// any object which meets the criterion is marked as -// PS_SOURCE_BRIGHTSTAR -# define MIN_DS 0.01 -bool mark_psf_source (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) -{ - int keep; - float dSX, dSY, SX, SY, SN; - float nSx, nSy, Chi; - - if (source->modelPSF == NULL) return (false); - - // if object has fitted peak above saturation, label as SATSTAR - // remember: fit does not use saturated pixels (masked) - if (source->modelPSF->params->data.F32[1] >= SATURATE) { - if (source->type == PS_SOURCE_PSFSTAR) { - psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n"); - } - source->type = PS_SOURCE_SATSTAR; - return (true); - } - if (source->type == PS_SOURCE_SATSTAR) { - psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); - source->type = PS_SOURCE_GOODSTAR; - } - - SN = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]; - SX = source->modelPSF->params->data.F32[4]; - SY = source->modelPSF->params->data.F32[5]; - dSX = source->modelPSF->dparams->data.F32[4]; - dSY = source->modelPSF->dparams->data.F32[5]; - Chi = source->modelPSF->chisq / source->modelPSF->nDOF; - - nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); - nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); - // sigma = 1 / SX - // dsx = 1 / (SX * SN) - // dsx_o = hypot (1/(SX*SN), MIN_DSX) - - // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] - keep = TRUE; - keep &= (fabs(nSx) < shapeNsigma); - keep &= (fabs(nSy) < shapeNsigma); - keep &= (SN > minSN); - keep &= (Chi < maxChi); - if (keep) { - if (source->type == PS_SOURCE_PSFSTAR) return (true); - source->type = PS_SOURCE_GOODSTAR; - return (true); - } - - if (source->type == PS_SOURCE_PSFSTAR) { - psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n"); - } - - // object appears to be small, suspected defect - if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { - source->type = PS_SOURCE_DEFECT; - return (false); - } - - // object appears to be large, suspected galaxy - if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { - source->type = PS_SOURCE_GALAXY; - return (false); - } - - // object appears to be extremely faint: what is this? - if (SN < minSN) { - source->type = PS_SOURCE_FAINTSTAR; - return (false); - } - - // these are pooly fitted, probable stars near other stars? - source->type = PS_SOURCE_POOR_FIT_PSF; - return (false); -} Index: unk/psphot/src/output.c =================================================================== --- /trunk/psphot/src/output.c (revision 4945) +++ (revision ) @@ -1,256 +1,0 @@ -# include "psphot.h" - -// output functions: we have several fixed modes (sx, obj, cmp) - -static int GetDophotType (psSource *source) { - - switch (source->type) { - - case PS_SOURCE_DEFECT: - case PS_SOURCE_SATURATED: - return (8); - - case PS_SOURCE_SATSTAR: - return (10); - - case PS_SOURCE_PSFSTAR: - case PS_SOURCE_GOODSTAR: - return (1); - - case PS_SOURCE_POOR_FIT_PSF: - return (7); - - case PS_SOURCE_FAIL_FIT_PSF: - case PS_SOURCE_FAINTSTAR: - return (4); - - case PS_SOURCE_GALAXY: - case PS_SOURCE_FAINT_GALAXY: - case PS_SOURCE_DROP_GALAXY: - case PS_SOURCE_FAIL_FIT_GAL: - case PS_SOURCE_POOR_FIT_GAL: - return (2); - - case PS_SOURCE_OTHER: - default: - return (0); - } - return (0); -} - -void output (psImageData *imdata, psMetadata *config, psArray *sources) { - - bool status; - - char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE"); - - if (outputMode == NULL) { - WriteSourcesText (imdata, config, sources); - exit (0); - } - - if (!strcasecmp (outputMode, "OBJ")) { - WriteSourcesOBJ (imdata, config, sources); - exit (0); - } - - if (!strcasecmp (outputMode, "SX")) { - WriteSourcesSX (imdata, config, sources); - exit (0); - } - - if (!strcasecmp (outputMode, "CMP")) { - WriteSourcesCMP (imdata, config, sources); - exit (0); - } - - if (!strcasecmp (outputMode, "CMF")) { - WriteSourcesCMF (imdata, config, sources); - exit (0); - } - - psAbort ("psphot", "unknown output mode %s", outputMode); -} - -bool WriteSourcesText (psImageData *imdata, psMetadata *config, psArray *sources) { - - bool status; - - char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); - - DumpImage (imdata->image, filename); - - // get these names from config? - DumpModelPSF (sources, "psfsources.dat"); - DumpModelFLT (sources, "fltsources.dat"); - DumpModelNULL (sources, "nullsources.dat"); - DumpMoments (sources, "moments.dat"); - - exit (0); -} - -// dophot-style output list with fixed line width -bool WriteSourcesOBJ (psImageData *imdata, psMetadata *config, psArray *sources) { - - int i, type; - bool status; - psModel *model; - psF32 *PAR, *dPAR; - float dmag, apMag, fitMag; - - char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); - 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 first - for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - model = source->modelPSF; - if (model == NULL) continue; - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); - dmag = dPAR[1] / PAR[1]; - - type = GetDophotType (source); - psLineInit (line); - psLineAdd (line, "%3d", type); - psLineAdd (line, "%8.2f", PAR[2]); - psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", 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", 32.0); - psLineAdd (line, "%8.3f", apMag); - psLineAdd (line, "%8.2f", apMag - fitMag); - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - return true; - -} - -// elixir/sextractor-style output list with fixed line width -bool WriteSourcesSX (psImageData *imdata, psMetadata *config, psArray *sources) { - - int i, type; - bool status; - psModel *model; - psF32 *PAR, *dPAR; - float dmag, apMag, fitMag; - - char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); - 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 first - for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelPSF; - if (model == NULL) continue; - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); - dmag = dPAR[1] / PAR[1]; - - // XXX EAM : fix this to match sextractor output - type = GetDophotType (source); - psLineInit (line); - psLineAdd (line, "%3d", type); - psLineAdd (line, "%8.2f", PAR[2]); - psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", 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", 32.0); - psLineAdd (line, "%8.3f", apMag); - psLineAdd (line, "%8.2f", apMag - fitMag); - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - return true; -} - -// elixir-style pseudo FITS table (header + ascii list) -bool WriteSourcesCMP (psImageData *imdata, psMetadata *config, psArray *sources) { - - int i, type; - bool status; - psModel *model; - psF32 *PAR, *dPAR; - float dmag, apMag, fitMag; - - char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); - 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 imdata->header to file - - // write sources with models first - for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelPSF; - if (model == NULL) continue; - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - - pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); - dmag = dPAR[1] / PAR[1]; - - // XXX EAM : fix this to match sextractor output - type = GetDophotType (source); - psLineInit (line); - psLineAdd (line, "%3d", type); - psLineAdd (line, "%8.2f", PAR[2]); - psLineAdd (line, "%8.2f", PAR[3]); - psLineAdd (line, "%8.3f", 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", 32.0); - psLineAdd (line, "%8.3f", apMag); - psLineAdd (line, "%8.2f", apMag - fitMag); - fwrite (line->line, 1, line->Nline, f); - } - fclose (f); - return true; -} - -// elixir-style FITS table output (header + table in 1st extension) -bool WriteSourcesCMF (psImageData *imdata, psMetadata *config, psArray *sources) { - - bool status; - - char *filename = psMetadataLookupPtr (&status, config, "OUTPUT"); - - // write imdata->header to file - - // write FITS table to file (use autocode tools) - - return true; -} Index: /trunk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 4946) +++ /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 4946) @@ -0,0 +1,43 @@ +# include "psphot.h" + +psArray *pmPeaksSigmaLimit (psImageData *imdata, psMetadata *config, psStats *sky) +{ + bool status = false; + float NSIGMA; + float SIGMA; + float threshold; + + // smooth the image + + psTimerStart ("psphot"); + + SIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_SIGMA"); + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA"); + + psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32); + psImageSmooth (smooth, SIGMA, NSIGMA); + psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot")); + + psTimerStart ("psphot"); + + // set peak threshold + NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); + threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold); + + // find the peaks in the smoothed image + psArray *peaks = pmFindImagePeaks (smooth, threshold); + if (peaks == NULL) psAbort ("find peaks", "no peaks found"); + psFree (smooth); + + // make this optional + // DumpPeaks (peaks, "peaks.dat"); + + psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); + + return (peaks); +} + +// In this function, we smooth the image, then search for the peaks +// Should we also subtract a super-binned image? (as an option?) +// We need to gracefully handle no source detections Index: /trunk/psphot/src/pmSourceUtils.c =================================================================== --- /trunk/psphot/src/pmSourceUtils.c (revision 4946) +++ /trunk/psphot/src/pmSourceUtils.c (revision 4946) @@ -0,0 +1,51 @@ +# include "psphot.h" + +bool pmSourceDefinePixels(psSource *mySource, + const psImageData *imdata, + psF32 x, + psF32 y, + psF32 Radius) +{ + psRegion *srcRegion; + + // Grab a subimage of the original image of size (2 * outerRadius). + srcRegion = psRegionSquare (x, y, Radius); + srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion); + + // use these when psImageSubset is updated + // psImage *subImage = psImageSubset(imdata->image, srcRegion); + // psImage *subImageMask = psImageSubset(imdata->mask, srcRegion); + // psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion); + + mySource->pixels = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); + mySource->noise = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); + mySource->mask = psImageSubset(imdata->mask, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1); + + psFree (srcRegion); + return(mySource); +} + +bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) { + + float obsSum = 0; + float fitSum = 0; + float sky = model->params->data.F32[0]; + + psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); + fitSum = modelFluxFunc (model->params); + + 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; + // fitSum += psModelEval (model, image, ix, iy) - sky; + } + } + if (obsSum <= 0) return false; + if (fitSum <= 0) return false; + + *fitMag = -2.5*log10(fitSum); + *obsMag = -2.5*log10(obsSum); + return (true); +} + Index: /trunk/psphot/src/psImageData.c =================================================================== --- /trunk/psphot/src/psImageData.c (revision 4946) +++ /trunk/psphot/src/psImageData.c (revision 4946) @@ -0,0 +1,24 @@ +# include "psphot.h" + +static void psImageDataFree (psImageData *imdata) { + + if (imdata == NULL) return; + + psFree (imdata->image); + psFree (imdata->header); + psFree (imdata->noise); + psFree (imdata->mask); + return; +} + +psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header) { + + psImageData *imdata = psAlloc(sizeof(psImageData)); + imdata->image = image; + imdata->header = header; + imdata->noise = noise; + imdata->mask = mask; + + p_psMemSetDeallocator(imdata, (psFreeFcn) psImageDataFree); + return (imdata); +} Index: /trunk/psphot/src/psLibUtils.c =================================================================== --- /trunk/psphot/src/psLibUtils.c (revision 4946) +++ /trunk/psphot/src/psLibUtils.c (revision 4946) @@ -0,0 +1,456 @@ +# include "psphot.h" + +// XXX EAM : these utility functions should be added back into PSLib + +static psHash *timers = NULL; + +bool psTimerClear (char *name) { + + bool status; + + if (name == NULL) return false; + + status = psHashRemove (timers, name); + return (status); +} + +void psTimerFree () { + + psFree (timers); + p_psTimeFinalize(); + return; +} + +// start/restart a named timer +bool psTimerStart (char *name) { + + psTime *start; + + if (timers == NULL) timers = psHashAlloc (16); + + start = psTimeGetTime (PS_TIME_UTC); + psHashAdd (timers, name, start); + psFree (start); + return (TRUE); +} + +// get current elapsed time on named timer +psF64 psTimerMark (char *name) { + + psTime *start; + psTime *mark; + psF64 delta; + + if (timers == NULL) return (0); + + start = psHashLookup (timers, name); + if (start == NULL) return (0); + + mark = psTimeGetTime (PS_TIME_UTC); + delta = psTimeDelta (mark, start); + psFree (mark); + // psFree (start); -- XXX is psHashLookup not psMemCopying? + + return (delta); +} + +// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) +// 2 = default, -v = 3, -vv = 4, -vvv = 5 +psS32 psLogArguments (int *argc, char **argv) { + + int N, level; + + // default log level is 2 + level = 2; + + // set in order, so that -vvv overrides -vv overrides -v + if ((N = psArgumentGet (*argc, argv, "-v"))) { + psArgumentRemove (N, argc, argv); + level = 3; + } + if ((N = psArgumentGet (*argc, argv, "-vv"))) { + psArgumentRemove (N, argc, argv); + level = 4; + } + if ((N = psArgumentGet (*argc, argv, "-vvv"))) { + psArgumentRemove (N, argc, argv); + level = 5; + } + + if ((N = psArgumentGet (*argc, argv, "-logfmt"))) { + if (*argc < N + 2) { + psAbort ("psLogArguments", "USAGE: -logfmt (format)"); + } + psArgumentRemove (N, argc, argv); + psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid + psArgumentRemove (N, argc, argv); + } + + // set the level, return the level + psLogSetLevel (level); + return (level); +} + +// set trace levels by facility +psS32 psTraceArguments (int *argc, char **argv) { + + int N; + + // argument format is: -trace (facil) (level) + while ((N = psArgumentGet (*argc, argv, "-trace"))) { + if (*argc < N + 3) { + psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); + } + psArgumentRemove (N, argc, argv); + psTraceSetLevel (argv[N], atoi(argv[N+1])); + psArgumentRemove (N, argc, argv); + psArgumentRemove (N, argc, argv); + } + if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) { + psTracePrintLevels (); + exit (2); + } + return (TRUE); +} + +// find the location of the specified argument +int psArgumentGet (int argc, char **argv, char *arg) { + + int i; + + for (i = 0; i < argc; i++) { + if (!strcmp(argv[i], arg)) + return (i); + } + + return ((int)NULL); +} + +// remove the specified argument (by location) +int psArgumentRemove (int N, int *argc, char **argv) { + + int i; + + if ((N != (int)NULL) && (N != 0)) { + (*argc)--; + for (i = N; i < *argc; i++) { + argv[i] = argv[i+1]; + } + } + + return (N); +} + +// alternate implementation of this function from pmObjects.c +psVector *psGetRowVectorFromImage(psImage *image, + psU32 row) +{ + PS_IMAGE_CHECK_NULL(image, NULL); + PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); + + psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); + memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32)); + return(tmpVector); +} + +// extract config informatin from config data or from image header, if specified +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_META_STR) { + 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_META_F32) { + value = item->data.F32; + // psFree (item); + return (value); + } + + *status = false; + psError(PS_ERR_UNKNOWN, true, "invalid item"); + return (0); +} + +void psImageSmooth (psImage *image, float sigma, float Nsigma) { + + int Nx, Ny, Npixel, Nrange; + float factor, g, s; + psVector *temp; + + // relevant terms + Nrange = sigma*Nsigma + 0.5; + Npixel = 2*Nrange + 1; + factor = -0.5/(sigma*sigma); + + Nx = image->numCols; + Ny = image->numRows; + + // generate gaussian + psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32); + for (int i = -Nrange; i < Nrange + 1; i++) { + gaussnorm->data.F32[i+Nrange] = exp (factor*i*i); + } + psF32 *gauss = &gaussnorm->data.F32[Nrange]; + + // smooth in X direction + temp = psVectorAlloc (Nx, PS_TYPE_F32); + for (int j = 0; j < Ny; j++) { + psF32 *vi = image->data.F32[j]; + psF32 *vo = temp->data.F32; + for (int i = 0; i < Nx; i++) { + g = s = 0; + for (int n = -Nrange; n < Nrange + 1; n++) { + if (i+n < 0) continue; + if (i+n >= Nx) continue; + s += gauss[n]*vi[i+n]; + g += gauss[n]; + } + vo[i] = s / g; + } + memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32)); + } + psFree (temp); + + // smooth in Y direction + temp = psVectorAlloc (image->numRows, PS_TYPE_F32); + for (int i = 0; i < Nx; i++) { + psF32 *vo = temp->data.F32; + psF32 **vi = image->data.F32; + for (int j = 0; j < Ny; j++) { + g = s = 0; + for (int n = -Nrange; n < Nrange + 1; n++) { + if (j+n < 0) continue; + if (j+n >= Ny) continue; + s += gauss[n]*vi[j+n][i]; + g += gauss[n]; + } + vo[j] = s / g; + } + // replace temp in image + for (int j = 0; j < Ny; j++) { + vi[j][i] = vo[j]; + } + } + psFree (temp); + psFree (gaussnorm); +} + +bool psImageInit (psImage *image,...) { + + va_list argp; + psU8 vU8; + psF32 vF32; + psF64 vF64; + + if (image == NULL) return (false); + + va_start (argp, image); + + switch (image->type.type) { + case PS_TYPE_U8: + vU8 = va_arg (argp, psU32); + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + image->data.U8[iy][ix] = vU8; + } + } + break; + + case PS_TYPE_F32: + vF32 = va_arg (argp, psF64); + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + image->data.F32[iy][ix] = vF32; + } + } + return (true); + + case PS_TYPE_F64: + vF64 = va_arg (argp, psF64); + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + image->data.F64[iy][ix] = vF64; + } + } + return (true); + + default: + psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type); + return (false); + } + return (false); +} + +// define a square region centered on the given coordinate +psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) { + psRegion *region; + region = psRegionAlloc (x - radius, x + radius + 1, + y - radius, y + radius + 1); + return (region); +} + +// set actual region based on image parameters: +// compensate for negative upper limits +// XXX this is inconsistent: the coordindates should always be in the parent +// frame, which means the negative values should subtract from Nx,Ny of +// the parent, not the child. but, we don't carry the dimensions of the +// parent in the psImage container. for now, us the child Nx,Ny +// force range to be on this subimage +// XXX EAM : this needs to be changes to use psRegion rather than psRegion* +psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) { + + // x0,y0, x1,y1 are in *parent* units + + if (out == NULL) { + out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1); + } else { + *out = *in; + } + // XXX these are probably wrong (see above) + if (out->x1 <= 0) { + out->x1 = image->col0 + image->numCols + out->x1; + } + if (out->y1 <= 0) { + out->y1 = image->row0 + image->numRows + out->y1; + } + + // force the lower-limits to be on the child + out->x0 = PS_MAX(image->col0, out->x0); + out->y0 = PS_MAX(image->row0, out->y0); + + // force the upper-limits to be on the child + out->x1 = PS_MIN(image->col0 + image->numCols, out->x1); + out->y1 = PS_MIN(image->row0 + image->numRows, out->y1); + return (out); +} + +// mask the area contained by the region +// the region is defined wrt the parent image +void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + if (ix + image->col0 < region->x0) continue; + if (ix + image->col0 >= region->x1) continue; + if (iy + image->row0 < region->y0) continue; + if (iy + image->row0 >= region->y1) continue; + if (logical_and) { + image->data.U8[iy][ix] &= maskValue; + } else { + image->data.U8[iy][ix] |= maskValue; + } + } + } +} + +// mask the area not contained by the region +// the region is defined wrt the parent image +void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + if (ix + image->col0 < region->x0) goto maskit; + if (ix + image->col0 >= region->x1) goto maskit; + if (iy + image->row0 < region->y0) goto maskit; + if (iy + image->row0 >= region->y1) goto maskit; + continue; + maskit: + if (logical_and) { + image->data.U8[iy][ix] &= maskValue; + } else { + image->data.U8[iy][ix] |= maskValue; + } + } + } +} + +// mask the area contained by the region +// the region is defined wrt the parent image +void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { + + double dx, dy, r2, R2; + + R2 = PS_SQR(radius); + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + dx = ix + image->col0 - x; + dy = iy + image->row0 - y; + r2 = PS_SQR(dx) + PS_SQR(dy); + if (r2 > R2) continue; + if (logical_and) { + image->data.U8[iy][ix] &= maskValue; + } else { + image->data.U8[iy][ix] |= maskValue; + } + } + } +} + +// mask the area contained by the region +// the region is defined wrt the parent image +void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { + + double dx, dy, r2, R2; + + R2 = PS_SQR(radius); + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + dx = ix + image->col0 - x; + dy = iy + image->row0 - y; + r2 = PS_SQR(dx) + PS_SQR(dy); + if (r2 < R2) continue; + if (logical_and) { + image->data.U8[iy][ix] &= maskValue; + } else { + image->data.U8[iy][ix] |= maskValue; + } + } + } +} + +psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { + + int nBin = (upper - lower) / delta; + + psVector *out = psVectorAlloc (nBin, type); + + for (int i = 0; i < nBin; i++) { + out->data.F64[i] = lower + i * delta; + } + + return (out); +} Index: /trunk/psphot/src/psLine.c =================================================================== --- /trunk/psphot/src/psLine.c (revision 4946) +++ /trunk/psphot/src/psLine.c (revision 4946) @@ -0,0 +1,46 @@ +# include "psphot.h" + +static void psLineFree (psLine *line) { + + if (line == NULL) return; + + psFree (line->line); + return; +} + +// allocate a psLine structrue +psLine *psLineAlloc (int Nline) { + + psLine *line; + line = psAlloc (sizeof(psLine)); + p_psMemSetDeallocator(line, (psFreeFcn) psLineFree); + + line->Nline = 0; + line->NLINE = Nline; + line->line = psAlloc (Nline); + return (line); +} + +bool psLineInit (psLine *line) { + if (line == NULL) return (false); + line->Nline = 0; + return (true); +} + +bool psLineAdd (psLine *line, char *format, ...) { + + int Nchar; + va_list ap; + + if (line == NULL) return (false); + + int nMax = line->NLINE - line->Nline; + + va_start (ap, format); + Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap); + line->Nline += PS_MIN (nMax - 1, Nchar); + + if (Nchar >= nMax) return (false); + return (true); +} + Index: /trunk/psphot/src/psMinimize.c =================================================================== --- /trunk/psphot/src/psMinimize.c (revision 4946) +++ /trunk/psphot/src/psMinimize.c (revision 4946) @@ -0,0 +1,250 @@ +# include "psphot.h" + +// XXX EAM : this file contains my alternate implementations of psMinimizeLMChi2 + +// XXX EAM this implementation of MinLM includes limits on params & dparams +psBool psMinimizeLMChi2_EAM(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func) +{ + PS_PTR_CHECK_NULL(min, NULL); + PS_VECTOR_CHECK_NULL(params, NULL); + PS_VECTOR_CHECK_EMPTY(params, NULL); + PS_PTR_CHECK_NULL(x, NULL); + PS_VECTOR_CHECK_NULL(y, NULL); + PS_VECTOR_CHECK_EMPTY(y, NULL); + PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL); + PS_PTR_CHECK_NULL(func, NULL); + + // this function has test and current values for several things + // the current best value is in lower case + // the next guess value is in upper case + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32); + psVector *dy = NULL; + psF64 Chisq = 0.0; + psF64 lambda = 0.001; + + psVector *beta_lim = NULL; + psVector *param_min = NULL; + psVector *param_max = NULL; + + // if we are provided a covar image, we expect to find these three vectors in first three rows + if (covar != NULL) { + beta_lim = psVectorAlloc (params->n, PS_TYPE_F32); + param_min = psVectorAlloc (params->n, PS_TYPE_F32); + param_max = psVectorAlloc (params->n, PS_TYPE_F32); + for (int i = 0; i < params->n; i++) { + beta_lim->data.F32[i] = covar->data.F64[0][i]; + param_min->data.F32[i] = covar->data.F64[1][i]; + param_max->data.F32[i] = covar->data.F64[2][i]; + } + psImageRecycle (covar, params->n, params->n, PS_TYPE_F64); + } + + + // why is this needed here??? the initial guess on params is provided by the user + Params = psVectorCopy (Params, params, PS_TYPE_F32); + + // the user provides the error or NULL. we need to convert + // to appropriate weights + dy = psVectorAlloc (y->n, PS_TYPE_F32); + if (yErr != NULL) { + for (int i = 0; i < dy->n; i++) { + if (yErr->data.F32[i] == 0.0) { + dy->data.F32[i] = 1.0; + // mask this? bad pixel, obviously... + } else { + dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]); + } + } + } else { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0; + } + } + + // calculate initial alpha and beta, set chisq (min->value) + min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); + if (isnan(min->value)) { + min->iter = min->maxIter; + return (false); + } + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { + p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), params, "params guess"); + } + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { + p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + // iterate until the tolerance is reached, or give up + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + + // set a new guess for Alpha, Beta, Params + p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda); + + // measure linear model prediction + psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda); + + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { + p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); + } + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { + p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + // calculate Chisq for new guess, update Alpha & Beta + Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func); + + // XXX EAM alternate convergence criterion: + // compare the delta (min->value - Chisq) with the + // expected delta from the linear model (dLinear) + // accept new guess (if improvement), or increase lambda + psF64 rho = (min->value - Chisq) / dLinear; + + psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho); + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { + p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + /* if (Chisq < min->value) { */ + if (rho > 0.0) { + min->lastDelta = (min->value - Chisq) / (dy->n - params->n); + min->value = Chisq; + alpha = psImageCopy (alpha, Alpha, PS_TYPE_F64); + beta = psVectorCopy (beta, Beta, PS_TYPE_F64); + params = psVectorCopy (params, Params, PS_TYPE_F32); + lambda *= 0.1; + } else { + lambda *= 10.0; + } + min->iter ++; + } + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); + + // construct & return the covariance matrix (if requested) + if (covar != NULL) { + p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0); + } + + // free the internal temporary data + psFree (alpha); + psFree (Alpha); + psFree (beta); + psFree (Beta); + psFree (Params); + psFree (dy); + + if (min->iter == min->maxIter) { + return (false); + } + return (true); +} + +// XXX EAM : can we use static copies of LUv, LUm, A? +psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, + psVector *Beta, + psVector *Params, + const psImage *alpha, + const psVector *beta, + const psVector *params, + const psVector *paramMask, + const psVector *beta_lim, + const psVector *params_min, + const psVector *params_max, + psF64 lambda) +{ + + # define USE_LU_DECOMP 1 + # if (USE_LU_DECOMP) + psVector *LUv = NULL; + psImage *LUm = NULL; + psImage *A = NULL; + psF32 det; + + // LU decomposition version + psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n"); + + // set new guess values (creates matrix A) + A = psImageCopy (NULL, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; + A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + // solve A*beta = Beta (Alpha = 1/A) + // these operations do not modify the input values (creates LUm, LUv) + LUm = psMatrixLUD (NULL, &LUv, A); + Beta = psMatrixLUSolve (Beta, LUm, beta, LUv); + Alpha = psMatrixInvert (Alpha, A, &det); + + # else + // gauss-jordan version + psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version"); + + // set new guess values (creates matrix A) + Beta = psVectorCopy (Beta, beta, PS_TYPE_F64); + Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); + for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; + Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); + } + + psGaussJordan (Alpha, Beta); + # endif + + // apply Beta to get new Params values + for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) continue; + // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; + // continue; + // compare Beta to beta limits + if (beta_lim != NULL) { + if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) { + Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]); + } + } + Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; + // compare new params to param limits + if (params_max != NULL) { + Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]); + } + if (params_min != NULL) { + Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]); + } + } + + # if (USE_LU_DECOMP) + psFree (A); + psFree (LUm); + psFree (LUv); + # endif + + return true; +} + Index: /trunk/psphot/src/psModulesUtils.c =================================================================== --- /trunk/psphot/src/psModulesUtils.c (revision 4946) +++ /trunk/psphot/src/psModulesUtils.c (revision 4946) @@ -0,0 +1,320 @@ +# include "psphot.h" + +// XXX EAM : these are my alternate implementations of psModule functions + +// this sets and clears bit 0x80 +bool pmSourceLocalSky_EAM (psSource *source, + psStatsOptions statsOptions, + psF32 Radius) +{ + + psImage *image = source->pixels; + psImage *mask = source->mask; + psPeak *peak = source->peak; + psRegion *srcRegion; + + srcRegion = psRegionSquare (peak->x, peak->y, Radius); + srcRegion = psRegionForImage (srcRegion, mask, srcRegion); + psImageMaskRegion (mask, srcRegion, OR, 0x80); + + psStats *myStats = psStatsAlloc(statsOptions); + myStats = psImageStats(myStats, image, mask, 0xff); + psImageMaskRegion (mask, srcRegion, AND, 0x7f); + + psF64 tmpF64; + p_psGetStatValue(myStats, &tmpF64); + + psFree (myStats); + psFree (srcRegion); + + if (isnan(tmpF64)) return (false); + source->moments = pmMomentsAlloc(); + source->moments->Sky = (psF32) tmpF64; + return (true); +} + +# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 +# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 +bool pmSourceFitModel_EAM(psSource *source, + psModel *model, + const bool PSF) +{ + PS_PTR_CHECK_NULL(source, false); + PS_PTR_CHECK_NULL(source->moments, false); + PS_PTR_CHECK_NULL(source->peak, false); + PS_PTR_CHECK_NULL(source->pixels, false); + PS_PTR_CHECK_NULL(source->mask, false); + PS_PTR_CHECK_NULL(source->noise, false); + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + psF32 Ro, ymodel; + + + // XXX EAM : is it necessary for the mask & noise to exist? the + // tests below could be conditions (!NULL) + + psModelFunc modelFunc = psModelFunc_GetFunction (model->type); + psModelLimits modelLimits = psModelLimits_GetFunction (model->type); + + psVector *params = model->params; + psVector *dparams = model->dparams; + psVector *paramMask = NULL; + + psVector *beta_lim = NULL; + psVector *params_min = NULL; + psVector *params_max = NULL; + + int nParams = PSF ? params->n - 4 : params->n; + psF32 So = params->data.F32[0]; + + // find the number of valid pixels + // XXX EAM : this loop and the loop below could just be one pass + // using the psArrayAdd and psVectorExtend functions + psS32 count = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + count++; + } + } + } + if (count < nParams + 1) { + psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); + return(false); + } + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(count); + psVector *y = psVectorAlloc(count, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(count, PS_TYPE_F32); + psS32 tmpCnt = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + if (source->mask->data.U8[i][j] == 0) { + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + // XXX: Convert i/j to image space: + // XXX EAM: coord order is (x,y) == (col,row) + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[tmpCnt] = (psPtr *) coord; + y->data.F32[tmpCnt] = source->pixels->data.F32[i][j]; + + ymodel = modelFunc (NULL, model->params, coord); + + // this test enhances the noise based on deviation from the model flux + Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So)); + yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro); + // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then + // the minimization function calculates sq() + tmpCnt++; + } + } + } + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + // PSF model only fits first 4 parameters, FLT model fits all + if (PSF) { + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + for (int i = 0; i < 4; i++) { + paramMask->data.U8[i] = 0; + } + for (int i = 4; i < paramMask->n; i++) { + paramMask->data.U8[i] = 1; + } + } + + psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); + modelLimits (&beta_lim, ¶ms_min, ¶ms_max); + for (int i = 0; i < params->n; i++) { + covar->data.F64[0][i] = beta_lim->data.F32[i]; + covar->data.F64[1][i] = params_min->data.F32[i]; + covar->data.F64[2][i] = params_max->data.F32[i]; + } + + psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if ((paramMask != NULL) && paramMask->data.U8[i]) continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // XXX EAM: we need to do something (give an error?) if rc is false + // XXX EAM: psMinimizeLMChi2 does not check convergence + + // XXX models can go insane: reject these + onPic &= (params->data.F32[2] >= source->pixels->col0); + onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (params->data.F32[3] >= source->pixels->row0); + onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + + // XXX EAM: save the resulting chisq, nDOF, nIter + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // XXX EAM get the Gauss-Newton distance for fixed model parameters + if (paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc); + for (int i = 0; i < dparams->n; i++) { + if (!paramMask->data.U8[i]) continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + } + + psFree(paramMask); + psFree(x); + psFree(y); + psFree(myMin); + + rc = (onPic && fitStatus); + return(rc); +} + +/****************************************************************************** +pmSourceMoments(source, radius): this function takes a subImage defined in the +psSource data structure, along with the peak location, and determines the +various moments associated with that peak. + +Requires the following to have been created: + psSource + psSource->peak + psSource->pixels + +XXX: The peak calculations are done in image coords, not subImage coords. + +XXX: mask values? +*****************************************************************************/ +# define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0) + +bool pmSourceMoments_EAM(psSource *source, + psF32 radius) +{ + PS_PTR_CHECK_NULL(source, NULL); + PS_PTR_CHECK_NULL(source->peak, NULL); + PS_PTR_CHECK_NULL(source->pixels, NULL); + PS_PTR_CHECK_NULL(source->noise, NULL); + PS_FLOAT_COMPARE(0.0, radius, NULL); + + // + // XXX: Verify the setting for sky if source->moments == NULL. + // + psF32 sky = 0.0; + if (source->moments == NULL) { + source->moments = pmMomentsAlloc(); + } else { + sky = source->moments->Sky; + } + + // + // Sum = SUM (z - sky) + // X1 = SUM (x - xc)*(z - sky) + // X2 = SUM (x - xc)^2 * (z - sky) + // XY = SUM (x - xc)*(y - yc)*(z - sky) + // + psF32 peakPixel = -PS_MAX_F32; + psS32 numPixels = 0; + psF32 Sum = 0.0; + psF32 X1 = 0.0; + psF32 Y1 = 0.0; + psF32 X2 = 0.0; + psF32 Y2 = 0.0; + psF32 XY = 0.0; + psF32 x = 0; + psF32 y = 0; + psF32 R2 = PS_SQR(radius); + + psF32 xPeak = source->peak->x; + psF32 yPeak = source->peak->y; + + // XXX why do I get different results for these two methods of finding Sx? + // XXX Sx, Sy would be better measured if we clip pixels close to sky + // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2? + // We loop through all pixels in this subimage (source->pixels), and for each + // pixel that is not masked, AND within the radius of the peak pixel, we + // proceed with the moments calculation. need to do two loops for a + // numerically stable result. first loop: get the sums. + // XXX EAM : mask == 0 is valid + + for (psS32 row = 0; row < source->pixels->numRows ; row++) { + for (psS32 col = 0; col < source->pixels->numCols ; col++) { + if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue; + + psF32 xDiff = col + source->pixels->col0 - xPeak; + psF32 yDiff = row + source->pixels->row0 - yPeak; + + if (!VALID_RADIUS(xDiff, yDiff)) continue; + + psF32 pDiff = source->pixels->data.F32[row][col] - sky; + + // XXX EAM : check for valid S/N in pixel + if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue; + + Sum += pDiff; + X1 += xDiff * pDiff; + Y1 += yDiff * pDiff; + XY += xDiff * yDiff * pDiff; + + X2 += PS_SQR(xDiff) * pDiff; + Y2 += PS_SQR(yDiff) * pDiff; + + peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); + numPixels++; + } + } + // XXX EAM - the limit is a bit arbitrary. make it user defined? + if ((numPixels < 3) || (Sum <= 0)) { + psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n"); + return (false); + } + + psTrace (".psModules.pmSourceMoments", 5, + "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", + sky, Sum, X1, Y1, X2, Y2, XY, numPixels); + + // + // first moment X = X1/Sum + xc + // second moment X = sqrt (X2/Sum - (X1/Sum)^2) + // Sxy = XY / Sum + // + x = X1/Sum; + y = Y1/Sum; + if ((fabs(x) > radius) || (fabs(y) > radius)) { + psTrace (".psModules.pmSourceMoments", 5, + "large centroid swing; invalid peak %d, %d\n", + source->peak->x, source->peak->y); + return (false); + } + + source->moments->x = x + xPeak; + source->moments->y = y + yPeak; + + source->moments->Sxy = XY/Sum - x*y; + source->moments->Sum = Sum; + source->moments->Peak = peakPixel; + source->moments->nPixels = numPixels; + + // XXX EAM : these values can be negative, so we need to limit the range + source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); + source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); + + psTrace (".psModules.pmSourceMoments", 4, + "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", + sky, Sum, source->moments->x, source->moments->y, + source->moments->Sx, source->moments->Sy, source->moments->Sxy); + + return(true); +} + +bool pmModelFitStatus (psModel *model) { + + bool status; + + psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type); + status = statusFunc (model); + + return (status); +} Index: unk/psphot/src/psUtils.c =================================================================== --- /trunk/psphot/src/psUtils.c (revision 4945) +++ (revision ) @@ -1,482 +1,0 @@ -# include "psphot.h" - -static psHash *timers = NULL; - -bool psTimerClear (char *name) { - - bool status; - - if (name == NULL) return false; - - status = psHashRemove (timers, name); - return (status); -} - -void psTimerFree () { - - psFree (timers); - p_psTimeFinalize(); - return; -} - -// start/restart a named timer -bool psTimerStart (char *name) { - - psTime *start; - - if (timers == NULL) timers = psHashAlloc (16); - - start = psTimeGetTime (PS_TIME_UTC); - psHashAdd (timers, name, start); - psFree (start); - return (TRUE); -} - -// get current elapsed time on named timer -psF64 psTimerMark (char *name) { - - psTime *start; - psTime *mark; - psF64 delta; - - if (timers == NULL) return (0); - - start = psHashLookup (timers, name); - if (start == NULL) return (0); - - mark = psTimeGetTime (PS_TIME_UTC); - delta = psTimeDelta (mark, start); - psFree (mark); - // psFree (start); -- XXX is psHashLookup not psMemCopying? - - return (delta); -} - -// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae) -// 2 = default, -v = 3, -vv = 4, -vvv = 5 -psS32 psLogArguments (int *argc, char **argv) { - - int N, level; - - // default log level is 2 - level = 2; - - // set in order, so that -vvv overrides -vv overrides -v - if ((N = get_argument (*argc, argv, "-v"))) { - remove_argument (N, argc, argv); - level = 3; - } - if ((N = get_argument (*argc, argv, "-vv"))) { - remove_argument (N, argc, argv); - level = 4; - } - if ((N = get_argument (*argc, argv, "-vvv"))) { - remove_argument (N, argc, argv); - level = 5; - } - - if ((N = get_argument (*argc, argv, "-logfmt"))) { - if (*argc < N + 2) { - psAbort ("psLogArguments", "USAGE: -logfmt (format)"); - } - remove_argument (N, argc, argv); - psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid - remove_argument (N, argc, argv); - } - - // set the level, return the level - psLogSetLevel (level); - return (level); -} - -// set trace levels by facility -psS32 psTraceArguments (int *argc, char **argv) { - - int N; - - // argument format is: -trace (facil) (level) - while ((N = get_argument (*argc, argv, "-trace"))) { - if (*argc < N + 3) { - psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)"); - } - remove_argument (N, argc, argv); - psTraceSetLevel (argv[N], atoi(argv[N+1])); - remove_argument (N, argc, argv); - remove_argument (N, argc, argv); - } - if ((N = get_argument (*argc, argv, "-trace-levels"))) { - psTracePrintLevels (); - exit (2); - } - return (TRUE); -} - -// find the location of the specified argument -int psArgumentGet (int argc, char **argv, char *arg) { - - int i; - - for (i = 0; i < argc; i++) { - if (!strcmp(argv[i], arg)) - return (i); - } - - return ((int)NULL); -} - -// remove the specified argument (by location) -int psArgumentRemove (int N, int *argc, char **argv) { - - int i; - - if ((N != (int)NULL) && (N != 0)) { - (*argc)--; - for (i = N; i < *argc; i++) { - argv[i] = argv[i+1]; - } - } - - return (N); -} - -// find the location of the specified argument -int get_argument (int argc, char **argv, char *arg) { - - int i; - - for (i = 0; i < argc; i++) { - if (!strcmp(argv[i], arg)) - return (i); - } - - return ((int)NULL); -} - -// remove the specified argument (by location) -int remove_argument (int N, int *argc, char **argv) { - - int i; - - if ((N != (int)NULL) && (N != 0)) { - (*argc)--; - for (i = N; i < *argc; i++) { - argv[i] = argv[i+1]; - } - } - - return (N); -} - -// alternate implementation of this function from pmObjects.c -psVector *psGetRowVectorFromImage(psImage *image, - psU32 row) -{ - PS_IMAGE_CHECK_NULL(image, NULL); - PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL); - - psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32); - memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32)); - return(tmpVector); -} - -// extract config informatin from config data or from image header, if specified -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_META_STR) { - 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_META_F32) { - value = item->data.F32; - // psFree (item); - return (value); - } - - *status = false; - psError(PS_ERR_UNKNOWN, true, "invalid item"); - return (0); -} - -void psImageSmooth (psImage *image, float sigma, float Nsigma) { - - int Nx, Ny, Npixel, Nrange; - float factor, g, s; - psVector *temp; - - // relevant terms - Nrange = sigma*Nsigma + 0.5; - Npixel = 2*Nrange + 1; - factor = -0.5/(sigma*sigma); - - Nx = image->numCols; - Ny = image->numRows; - - // generate gaussian - psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32); - for (int i = -Nrange; i < Nrange + 1; i++) { - gaussnorm->data.F32[i+Nrange] = exp (factor*i*i); - } - psF32 *gauss = &gaussnorm->data.F32[Nrange]; - - // smooth in X direction - temp = psVectorAlloc (Nx, PS_TYPE_F32); - for (int j = 0; j < Ny; j++) { - psF32 *vi = image->data.F32[j]; - psF32 *vo = temp->data.F32; - for (int i = 0; i < Nx; i++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (i+n < 0) continue; - if (i+n >= Nx) continue; - s += gauss[n]*vi[i+n]; - g += gauss[n]; - } - vo[i] = s / g; - } - memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32)); - } - psFree (temp); - - // smooth in Y direction - temp = psVectorAlloc (image->numRows, PS_TYPE_F32); - for (int i = 0; i < Nx; i++) { - psF32 *vo = temp->data.F32; - psF32 **vi = image->data.F32; - for (int j = 0; j < Ny; j++) { - g = s = 0; - for (int n = -Nrange; n < Nrange + 1; n++) { - if (j+n < 0) continue; - if (j+n >= Ny) continue; - s += gauss[n]*vi[j+n][i]; - g += gauss[n]; - } - vo[j] = s / g; - } - // replace temp in image - for (int j = 0; j < Ny; j++) { - vi[j][i] = vo[j]; - } - } - psFree (temp); - psFree (gaussnorm); -} - -bool psImageInit (psImage *image,...) { - - va_list argp; - psU8 vU8; - psF32 vF32; - psF64 vF64; - - if (image == NULL) return (false); - - va_start (argp, image); - - switch (image->type.type) { - case PS_TYPE_U8: - vU8 = va_arg (argp, psU32); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.U8[iy][ix] = vU8; - } - } - break; - - case PS_TYPE_F32: - vF32 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F32[iy][ix] = vF32; - } - } - return (true); - - case PS_TYPE_F64: - vF64 = va_arg (argp, psF64); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - image->data.F64[iy][ix] = vF64; - } - } - return (true); - - default: - psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type); - return (false); - } - return (false); -} - -// define a square region centered on the given coordinate -psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) { - psRegion *region; - region = psRegionAlloc (x - radius, x + radius + 1, - y - radius, y + radius + 1); - return (region); -} - -// set actual region based on image parameters: -// compensate for negative upper limits -// XXX this is inconsistent: the coordindates should always be in the parent -// frame, which means the negative values should subtract from Nx,Ny of -// the parent, not the child. but, we don't carry the dimensions of the -// parent in the psImage container. for now, us the child Nx,Ny -// force range to be on this subimage -// XXX EAM : this needs to be changes to use psRegion rather than psRegion* -psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) { - - // x0,y0, x1,y1 are in *parent* units - - if (out == NULL) { - out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1); - } else { - *out = *in; - } - // XXX these are probably wrong (see above) - if (out->x1 <= 0) { - out->x1 = image->col0 + image->numCols + out->x1; - } - if (out->y1 <= 0) { - out->y1 = image->row0 + image->numRows + out->y1; - } - - // force the lower-limits to be on the child - out->x0 = PS_MAX(image->col0, out->x0); - out->y0 = PS_MAX(image->row0, out->y0); - - // force the upper-limits to be on the child - out->x1 = PS_MIN(image->col0 + image->numCols, out->x1); - out->y1 = PS_MIN(image->row0 + image->numRows, out->y1); - return (out); -} - -// mask the area contained by the region -// the region is defined wrt the parent image -void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) continue; - if (ix + image->col0 >= region->x1) continue; - if (iy + image->row0 < region->y0) continue; - if (iy + image->row0 >= region->y1) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} - -// mask the area not contained by the region -// the region is defined wrt the parent image -void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) { - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (ix + image->col0 < region->x0) goto maskit; - if (ix + image->col0 >= region->x1) goto maskit; - if (iy + image->row0 < region->y0) goto maskit; - if (iy + image->row0 >= region->y1) goto maskit; - continue; - maskit: - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} - -// mask the area contained by the region -// the region is defined wrt the parent image -void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 > R2) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} - -// mask the area contained by the region -// the region is defined wrt the parent image -void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) { - - double dx, dy, r2, R2; - - R2 = PS_SQR(radius); - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - dx = ix + image->col0 - x; - dy = iy + image->row0 - y; - r2 = PS_SQR(dx) + PS_SQR(dy); - if (r2 < R2) continue; - if (logical_and) { - image->data.U8[iy][ix] &= maskValue; - } else { - image->data.U8[iy][ix] |= maskValue; - } - } - } -} - -psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) { - - int nBin = (upper - lower) / delta; - - psVector *out = psVectorAlloc (nBin, type); - - for (int i = 0; i < nBin; i++) { - out->data.F64[i] = lower + i * delta; - } - - return (out); -} Index: unk/psphot/src/psphot-utils.c =================================================================== --- /trunk/psphot/src/psphot-utils.c (revision 4945) +++ (revision ) @@ -1,464 +1,0 @@ -# include "psphot.h" - -// dump the peaks to an output file -bool DumpPeaks (psArray *sources, char *filename) -{ - - int i; - FILE *f; - psPeak *peak; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpPeaks", 3, "can't open output file for peaks%s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - peak = (psPeak *) sources->data[i]; - if (peak == NULL) continue; - fprintf (f, "%5d %5d %7.1f\n", - peak->x, peak->y, peak->counts); - } - fclose (f); - return true; -} - -// dump the moments to an output file -bool DumpMoments (psArray *sources, char *filename) -{ - - int i; - FILE *f; - psSource *source; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpMoments", 3, "can't open output file for moments%s\n", filename); - return false; - } - - for (i = 0; i < sources->n; i++) { - source = (psSource *) sources->data[i]; - if (source->moments == NULL) continue; - 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); - return true; -} - -bool DumpPSFTestData (pmPSF_Test *test, char *filename) -{ - - FILE *f = fopen (filename, "w"); - - for (int j = 0; j < test->sources->n; j++) { - - psSource *source = test->sources->data[j]; - psModel *model = test->modelPSF->data[j]; - - float sky = 0; - if (model != NULL) sky = model->params->data.F32[0]; - - fprintf (f, "%3d %6.1f %6.1f %7.4f %7.4f %8.4f\n", - j, source->moments->x, source->moments->y, - test->metric->data.F64[j], test->fitMag->data.F64[j], sky); - } - fclose (f); - return true; -} - -// dump the sources to an output file -bool DumpModelPSF (psArray *sources, char *filename) -{ - - double dP, flux; - int i, j; - FILE *f; - psVector *params; - psVector *dparams; - psModel *model; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpObjects", 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++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelPSF; - if (model == NULL) continue; - - // valid source types for this function - if (source->type == PS_SOURCE_SATSTAR) goto valid; - if (source->type == PS_SOURCE_PSFSTAR) goto valid; - if (source->type == PS_SOURCE_GOODSTAR) goto valid; - continue; - - valid: - params = model->params; - dparams = model->dparams; - - dP = 0; - dP += PS_SQR(dparams[0].data.F32[2]); - dP += PS_SQR(dparams[0].data.F32[3]); - dP = sqrt (dP); - - psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); - flux = modelFluxFunc (params); - - fprintf (f, "%7.1f %7.1f %5.1f %7.3f %7.4f %7.4f ", - params[0].data.F32[2], - params[0].data.F32[3], - params[0].data.F32[0], - -2.5*log10(flux), - (dparams[0].data.F32[1]/params[0].data.F32[1]), - dP); - for (j = 0; j < model->params->n - 4; j++) { - fprintf (f, "%9.6f ", params[0].data.F32[j+4]); - } - fprintf (f, " : "); - for (j = 0; j < model->params->n - 4; j++) { - fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]); - } - fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", - source[0].type, - 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 DumpModelFLT (psArray *sources, char *filename) -{ - - double dP; - int i, j; - FILE *f; - psVector *params; - psVector *dparams; - psModel *model; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpObjects", 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++) { - psSource *source = (psSource *) sources->data[i]; - model = (psModel *) source->modelFLT; - if (model == NULL) continue; - if (source->type == PS_SOURCE_GALAXY) goto valid; - continue; - - valid: - params = model->params; - dparams = model->dparams; - - // XXX these are hardwired for SGAUSS : this should be pushed into the - // model functions as an abstract function - dP = 0; - dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]); - dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]); - dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]); - dP = sqrt (dP); - - fprintf (f, "%7.1f %7.1f %5.1f %7.3f %7.4f %7.4f ", - params[0].data.F32[2], - params[0].data.F32[3], - params[0].data.F32[0], - -2.5*log10(params[0].data.F32[1]), - (dparams[0].data.F32[1]/params[0].data.F32[1]), - dP); - - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", params[0].data.F32[j]); - } - fprintf (f, " : "); - for (j = 4; j < model->params->n; j++) { - fprintf (f, "%9.6f ", dparams[0].data.F32[j]); - } - fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", - source[0].type, 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 DumpModelNULL (psArray *sources, char *filename) -{ - - int i; - FILE *f; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); - return false; - } - - psMoments *empty = pmMomentsAlloc (); - - // write sources with models first - for (i = 0; i < sources->n; i++) { - psSource *source = (psSource *) sources->data[i]; - - // skip these sources (in PSF or FLT) - if (source->type == PS_SOURCE_GALAXY) continue; - if (source->type == PS_SOURCE_PSFSTAR) continue; - if (source->type == PS_SOURCE_SATSTAR) continue; - if (source->type == PS_SOURCE_GOODSTAR) continue; - - if (source->moments == NULL) { - 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, - empty->x, empty->y, - empty->Sx, empty->Sy, - empty->Sum, empty->Peak, - empty->Sky, empty->SN, - empty->nPixels, source->type); - } else { - 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); - return true; -} - -// dump the models to an output file -bool DumpModelFits (psArray *models, char *filename) -{ - - int i, j; - FILE *f; - psVector *params; - psVector *dparams; - psModel *model; - - f = fopen (filename, "w"); - if (f == NULL) { - psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); - return false; - } - - // this loop is valid for GAUSS, PGAUSS: - for (i = 0; i < models->n; i++) { - model = (psModel *) models->data[i]; - if (model == NULL) continue; - params = model->params; - dparams = model->dparams; - fprintf (f, "%7.1f %7.1f %5.1f %7.1f ", - params[0].data.F32[2], params[0].data.F32[3], - params[0].data.F32[0], params[0].data.F32[1]); - for (j = 0; j < model->params->n - 4; j++) { - fprintf (f, "%9.6f ", params[0].data.F32[j+4]); - } - for (j = 0; j < model->params->n - 4; j++) { - fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]); - } - fprintf (f, ": %7.1f (%d %d)\n", model[0].chisq, model[0].nDOF, model[0].nIter); - } - fclose (f); - return true; -} - -bool DumpImage (psImage *image, char *filename) { - - unlink (filename); - psFits *fits = psFitsAlloc (filename); - psFitsWriteImage (fits, NULL, image, 0, NULL); - psFree (fits); - return true; -} - -bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) { - - float obsSum = 0; - float fitSum = 0; - float sky = model->params->data.F32[0]; - - psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); - fitSum = modelFluxFunc (model->params); - - 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; - // fitSum += psModelEval (model, image, ix, iy) - sky; - } - } - if (obsSum <= 0) return false; - if (fitSum <= 0) return false; - - *fitMag = -2.5*log10(fitSum); - *obsMag = -2.5*log10(obsSum); - return (true); -} - -# if (0) -void pmSourceMaskSaturated (psSource *source, float saturate) { - - psImage *image = source->pixels; - psImage *mask = source->mask; - - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - if (image->data.F32[iy][ix] < saturate) continue; - mask->data.U8[iy][ix] |= 0x02; - } - } -} -# endif - -// create a subset of sources for the PS_SOURCE_PSFSTAR entries -// XXX deprecated -psArray *SelectPSFStars (psArray *sources) { - - psArray *stars; - psSource *source; - int i; - - stars = psArrayAlloc (sources->n); - stars->n = 0; - - for (i = 0; i < sources->n; i++) { - source = (psSource *) sources->data[i]; - if (source->type != PS_SOURCE_PSFSTAR) continue; - psArrayAdd (stars, 200, source); - } - - return (stars); -} - -/* this code is in pmObjects.h / psModules:psEllipse */ -# if (0) -typedef struct { - double major; - double minor; - double theta; -} EllipseAxes; - -typedef struct { - double x2; - double y2; - double xy; -} EllipseMoments; - -typedef struct { - double sx; - double sy; - double sxy; -} EllipseShape; - -EllipseAxes EllipseMomentsToAxes (EllipseMoments moments) { - - EllipseAxes axes; - - double f = sqrt (0.25*PS_SQR(moments.x2 - moments.y2) + PS_SQR(moments.sy)); - - axes.major = sqrt (0.5*(moments.x2 + moments.y2) + f); - axes.minor = sqrt (0.5*(moments.x2 + moments.y2) - f) - axes.theta = atan2 (2*moments.xy, moments.x2 - moments.y2) / 2; - // theta in radians - - return (axes); -} - -EllipseShape EllipseMomentsToShape (EllipseMoments moments) { - - EllipseShape shape; - - double f = sqrt(0.25*PS_SQR(axes.x2 - axes.y2) + PS_SQR(axes.sxy)); - - shape.sx = moments.y2 / f; - shape.sy = moments.x2 / f; - shape.sxy = -2.0*moments.xy / f; - // XXX - these are probably wrong - - return (shape); -} - -EllipseMoments EllipseShapeToMoments (EllipseShape shape) { - - EllipseMoments moments; - - return (moments); -} - -EllipseAxes EllipseShapeToAxes (EllipseShape shape) { - - EllipseAxes axes; - - - return (axes); -} - -EllipseShape EllipseAxesToShape (EllipseAxes axes) { - - EllipseShape shape; - - shape.sx = PS_SQR(cos(axes.theta) / axes.major) + PS_SQR(sin(axes.theta) / axes.minor); - shape.sy = PS_SQR(sin(axes.theta) / axes.major) + PS_SQR(cos(axes.theta) / axes.minor); - shape.sxy = cos(axes.theta)*sin(axes.theta)*(1.0/PS_SQR(axes.major) - 1.0/PS_SQR(axes.minor)); - - return (shape); -} - -EllipseMoments EllipseAxesToMoments (EllipseAxes axes) { - - EllipseMoments moments; - - return (moments); -} - -# endif - -static void psImageDataFree (psImageData *imdata) { - - if (imdata == NULL) return; - - psFree (imdata->image); - psFree (imdata->header); - psFree (imdata->noise); - psFree (imdata->mask); - return; -} - -psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header) { - - psImageData *imdata = psAlloc(sizeof(psImageData)); - imdata->image = image; - imdata->header = header; - imdata->noise = noise; - imdata->mask = mask; - - p_psMemSetDeallocator(imdata, (psFreeFcn) psImageDataFree); - return (imdata); -} - Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 4945) +++ /trunk/psphot/src/psphot.c (revision 4946) @@ -9,39 +9,43 @@ pmPSF *psf = NULL; psStats *sky = NULL; + pmPSFClump psfClump; - config = load_args (&argc, argv); + config = psphotArguments (&argc, argv); // load input data (image and config) // create or load mask and noise images // we have memory leaks here -- may be from psMetadata - imdata = setup (config); + imdata = psphotSetup (config); // measure image stats for initial guess - sky = image_stats (imdata, config); + sky = psphotImageStats (imdata, config); // find the peaks in the image - peaks = find_peaks (imdata, config, sky); + peaks = pmPeaksSigmaLimit (imdata, config, sky); // construct sources and measure basic stats - sources = source_moments (imdata, config, peaks); + sources = psphotSourceStats (imdata, config, peaks); // classify sources based on moments, brightness - basic_classes (sources, config); + // XXX - this has 1 leak + psfClump = pmSourcePSFClump (sources, config); + + // group into STAR, COSMIC, GALAXY, SATURATED + pmSourceRoughClass (sources, config, psfClump); // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, by_SN); + sources = psArraySort (sources, psphotSortBySN); // use bright stellar objects to measure PSF - psf = choose_psf_model (config, sources, sky); + psf = psphotChoosePSF (config, sources, sky); // test PSF on all except SATURATE and DEFECT (artifacts) - apply_psf_model (imdata, config, sources, psf, sky); + psphotApplyPSF (imdata, config, sources, psf, sky); // fit extended objects with galaxy models - fit_galaxies (imdata, config, sources, sky); + psphotFitGalaxies (imdata, config, sources, sky); // write out data in appropriate format - output (imdata, config, sources); - + psphotOutput (imdata, config, sources); exit (0); } Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 4945) +++ /trunk/psphot/src/psphot.h (revision 4946) @@ -38,4 +38,10 @@ } psLine; +enum { + PSPHOT_MASK_INVALID = 0x01, + PSPHOT_MASK_SATURATED = 0x02, + PSPHOT_MASK_KEEP = 0x08, +} MaskValues; + // used by mask operations # define AND true @@ -45,119 +51,108 @@ // top-level psphot functions -psImageData *setup (psMetadata *config); -psStats *image_stats (psImageData *imdata, psMetadata *config); -psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky); -psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks); -pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *sky); -bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); -bool test_psf_scatter (psArray *sources); -bool mark_psf_sources (psArray *sources, psMetadata *config); -bool subtract_psf_sources (psArray *sources); -bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats); -bool subtract_galaxies (psArray *sources, psMetadata *config); -bool subtract_psf_source (psSource *source); -bool mark_psf_source (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); -bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf); -bool basic_classes (psArray *sources, psMetadata *config); +psMetadata *psphotArguments (int *argc, char **argv); +psImageData *psphotSetup (psMetadata *config); +psStats *psphotImageStats (psImageData *imdata, psMetadata *config); +psArray *pmPeaksSigmaLimit (psImageData *imdata, psMetadata *config, psStats *sky); +psArray *psphotSourceStats (psImageData *imdata, psMetadata *config, psArray *allpeaks); +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky); +bool psphotApplyPSF (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +bool psphotFitGalaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats); +void psphotOutput (psImageData *imdata, psMetadata *config, psArray *sources); + +bool psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE); +bool psphotSubtractPSF (psSource *source); +int psphotSortBySN (const void **a, const void **b); + +// psImageData functions +psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); // psf utilities -pmPSF *pmPSFAlloc (psModelType type); -pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName); -pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius); -bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask); -psModel *psModelFromPSF (psModel *model, pmPSF *psf); -bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask); -bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS); +pmPSF *pmPSFAlloc (psModelType type); +pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName); +pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius); +bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask); +psModel *psModelFromPSF (psModel *model, pmPSF *psf); +bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask); +bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS); -// output options -bool DumpPeaks (psArray *sources, char *filename); -bool DumpMoments (psArray *sources, char *filename); -bool DumpModelFits (psArray *models, char *filename); -bool DumpModelPSF (psArray *sources, char *filename); -bool DumpModelFLT (psArray *sources, char *filename); -bool DumpModelNULL (psArray *sources, char *filename); -bool DumpImage (psImage *image, char *filename); -bool DumpPSFTestData (pmPSF_Test *test, char *filename); -void output (psImageData *imdata, psMetadata *config, psArray *sources); -bool WriteSourcesText (psImageData *imdata, psMetadata *config, psArray *sources); -bool WriteSourcesOBJ (psImageData *imdata, psMetadata *config, psArray *sources); -bool WriteSourcesSX (psImageData *imdata, psMetadata *config, psArray *sources); -bool WriteSourcesCMP (psImageData *imdata, psMetadata *config, psArray *sources); -bool WriteSourcesCMF (psImageData *imdata, psMetadata *config, psArray *sources); -bool init_line (char *line, int *Nline); -bool add_to_line (char *line, int *Nline, char *format, ...); +// output functions +bool pmSourcesWriteText (psImageData *imdata, char *filename, psArray *sources); +bool pmSourcesWriteOBJ (psImageData *imdata, char *filename, psArray *sources); +bool pmSourcesWriteCMP (psImageData *imdata, char *filename, psArray *sources); +bool pmSourcesWriteCMF (psImageData *imdata, char *filename, psArray *sources); +bool pmSourcesWriteSX (psImageData *imdata, char *filename, psArray *sources); +int pmSourcesDophotType (psSource *source); +bool psPeaksWriteText (psArray *sources, char *filename); +bool psMomentsWriteText (psArray *sources, char *filename); +bool psModelWritePSFs (psArray *sources, char *filename); +bool psModelWriteFLTs (psArray *sources, char *filename); +bool psModelWriteNULLs (psArray *sources, char *filename); -// psphot utilities -int usage (); -psMetadata *load_args (int *argc, char **argv); -int by_SN (const void **a, const void **b); -psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid); -psArray *SelectPSFStars (psArray *sources); -bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius); -bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius); -bool pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF); -bool pmSourceMoments_EAM(psSource *source, psF32 radius); -bool pmModelFitStatus (psModel *model); -psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, psVector *Beta, psVector *Params, const psImage *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda); -psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func); -psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda); -psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header); +// psModule extra utilities +bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius); +bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius); +bool pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF); +bool pmSourceMoments_EAM(psSource *source, psF32 radius); +bool pmModelFitStatus (psModel *model); +int pmSourceDophotType (psSource *source); -// fitsource utilities -bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y); -psMetadata *fs_args (int *argc, char **argv); -int fs_usage (); +// minimize +psBool p_psMinLM_GuessABP_EAM (psImage *Alpha, psVector *Beta, psVector *Params, const psImage *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda); +psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func); +psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda); -// image mask functions -void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); -void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); -void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); +// psLib extra utilities +bool psTimerStart (char *name); +void psTimerFree (); +bool psTimerClear (char *name); +psF64 psTimerMark (char *name); +psS32 psLogArguments (int *argc, char **argv); +psS32 psTraceArguments (int *argc, char **argv); +int psArgumentGet (int argc, char **argv, char *arg); +int psArgumentRemove (int N, int *argc, char **argv); +psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); +psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); +void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); +void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue); +void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); +void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue); +psVector *psGetRowVectorFromImage(psImage *image, psU32 row); // basic image functions -bool psImageInit (psImage *image,...); -void psImageSmooth (psImage *image, float sigma, float Nsigma); -psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in); -psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius); - -// misc global functions -psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row); -psVector *psGetRowVectorFromImage(psImage *image, psU32 row); - -// utilities -bool psTimerStart (char *name); -void psTimerFree (); -bool psTimerClear (char *name); -psF64 psTimerMark (char *name); -psS32 psLogArguments (int *argc, char **argv); -psS32 psTraceArguments (int *argc, char **argv); -int get_argument (int argc, char **argv, char *arg); -int remove_argument (int N, int *argc, char **argv); -psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name); -psVector *psVectorCreate (double lower, double upper, double delta, psElemType type); +bool psImageInit (psImage *image,...); +void psImageSmooth (psImage *image, float sigma, float Nsigma); +psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in); +psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius); +int psphotSaveImage (psMetadata *header, psImage *image, char *filename); // psLine functions -psLine *psLineAlloc (int Nline); -bool psLineInit (psLine *line); -bool psLineAdd (psLine *line, char *format, ...); +psLine *psLineAlloc (int Nline); +bool psLineInit (psLine *line); +bool psLineAdd (psLine *line, char *format, ...); // polynomial functions -psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y); -psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm); -psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); -psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type); -void psPolynomial2DDump (psPolynomial2D *poly); -psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); -psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); -psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y); +psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y); +psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm); +psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); +psPolynomial2D *Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type); +void psPolynomial2DDump (psPolynomial2D *poly); +psPolynomial2D *RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); +psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr); +psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y); +psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm); +void psPolynomial1DDump (psPolynomial1D *poly); +psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly); +psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x); +psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type); +psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr); -psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm); -void psPolynomial1DDump (psPolynomial1D *poly); -psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly); -psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x); -psPolynomial1D* Polynomial1DAlloc(psS32 nOrder, psPolynomialType type); -psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr); +// fitsource +bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y); +psMetadata *fs_args (int *argc, char **argv); +int fs_usage (); +bool DumpImage (psImage *image, char *filename); // deprecated -void pmSourceMaskRegion (psSource *source, psRegion *region); -void pmSourceMaskSaturated (psSource *source, float saturate); +// void pmSourceMaskRegion (psSource *source, psRegion *region); +//psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row); Index: /trunk/psphot/src/psphotApplyPSF.c =================================================================== --- /trunk/psphot/src/psphotApplyPSF.c (revision 4946) +++ /trunk/psphot/src/psphotApplyPSF.c (revision 4946) @@ -0,0 +1,93 @@ +# include "psphot.h" + +// fit psf model to all objects +// PSFSTAR objects will be refitted +// run this function to a specific flux limit? + +bool psphotApplyPSF (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x; + float y; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + // we may set this differently here from the value used to mark likely saturated stars + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + + float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + // set the object surface-brightness limit for fitted pixels + float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // this function specifies the radius at this the model hits the given flux + psModelRadius modelRadius = psModelRadius_GetFunction (psf->type); + + for (int i = 0; i < sources->n; i++) { + + psSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (source->type == PS_SOURCE_DEFECT) continue; + if (source->type == PS_SOURCE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + psModel *modelFLT = pmSourceModelGuess (source, psf->type); + + // set PSF parameters for this model + psModel *model = psModelFromPSF (modelFLT, psf); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + psFree (modelFLT); + + // set the fit radius based on the object flux limit and the model + // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) + model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) { + psAbort ("apply_psf_model", "error in radius"); + } + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > OUTER_RADIUS) { + // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, imdata, x, y, model->radius); + } + + // fit PSF model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP); + status = pmSourceFitModel (source, model, true); + psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP); + if (!status || (model->params->data.F32[1] < 0)) { + psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter); + source->type = PS_SOURCE_FAIL_FIT_PSF; // better choice? + psFree (model); + continue; + } + + // model succeeded : tentatively keep it + source->modelPSF = model; + Niter += model[0].nIter; + Nfit ++; + + // is it a good model? + psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION); + if (psphotSubtractPSF (source)) { + Nsub ++; + } + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub); + return (true); +} Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 4946) +++ /trunk/psphot/src/psphotArguments.c (revision 4946) @@ -0,0 +1,73 @@ +# include "psphot.h" + +static void usage (void) ; + +psMetadata *psphotArguments (int *argc, char **argv) { + + int N, Nfail; + int mode = PS_META_STR | PS_META_REPLACE; + + // basic pslib options + fprintf (stderr, "starting...\n"); + psLogSetFormat ("M"); + psLogArguments (argc, argv); + psTraceArguments (argc, argv); + + // optional mask image - add to config + char *mask = NULL; + if ((N = psArgumentGet (*argc, argv, "-mask"))) { + psArgumentRemove (N, argc, argv); + mask = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + // optional noise image - add to config + char *noise = NULL; + if ((N = psArgumentGet (*argc, argv, "-noise"))) { + psArgumentRemove (N, argc, argv); + noise = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + // optional output residual image - add to config + char *resid = NULL; + if ((N = psArgumentGet (*argc, argv, "-resid"))) { + psArgumentRemove (N, argc, argv); + resid = psStringCopy (argv[N]); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 4) usage (); + + // load config information + psMetadata *config = psMetadataAlloc (); + psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_META_MULTI, "folder for psf model entries", NULL); + config = psMetadataParseConfig (config, &Nfail, argv[3], FALSE); + fprintf (stderr, "loaded config...\n"); + + // identify input image & optional noise & mask images + // command-line entries override config-file entries + psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); + psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]); + + if (mask != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask); + } + if (noise != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "NOISE_IMAGE", mode, "", noise); + } + if (resid != NULL) { + psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid); + } + return (config); +} + +static void usage (void) { + + fprintf (stderr, "USAGE: psphot (image) (output) (config)\n"); + fprintf (stderr, "options: \n"); + fprintf (stderr, " -mask (filename)\n"); + fprintf (stderr, " -noise (filename)\n"); + fprintf (stderr, " -resid (filename)\n"); + exit (2); +} Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 4946) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 4946) @@ -0,0 +1,94 @@ +# include "psphot.h" + +// test PSF models and select best option + +pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) +{ + bool status; + char *modelName; + pmPSF_Test *test = NULL; + psArray *stars = NULL; + psMetadataItem *item = NULL; + + // array to store candidate PSF stars + stars = psArrayAlloc (sources->n); + stars->n = 0; + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; i < sources->n; i++) { + psSource *source = sources->data[i]; + if (source->type != PS_SOURCE_PSFSTAR) continue; + psArrayAdd (stars, 200, source); + } + psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n); + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + + // get the list pointers for the PSF_MODEL entries + psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + + psList *list = (psList *) mdi->data.list; + psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); + + // set up an array to store the test results + psArray *tests = psArrayAlloc (list->size); + // XXX EAM : new value : psArray *tests = psArrayAlloc (list->n); + + // test each model option listed in config + for (int i = 0; i < tests->n; i++) { + + item = psListGetAndIncrement (iter); + modelName = item->data.V; + + tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS); + psFree (modelName); + psFree (item); + } + psFree (iter); + // psFree (list); XXX EAM - is list freed with iter? + psFree (stars); + + // select the best of the test models + // here we are using the clippedStdev on the metric as the indicator + test = tests->data[0]; + int bestN = 0; + float bestM = test->dApResid; + for (int i = 1; i < tests->n; i++) { + test = tests->data[i]; + float M = test->dApResid; + if (M < bestM) { + bestM = M; + bestN = i; + } + } + + // keep only the selected test: + test = tests->data[bestN]; + modelName = psModelGetType (test->modelType); + psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid); + + // XXX I am having trouble tracking down the double frees vs leaks + # if (0) + // we are keeping the modelPSF fits from the PSF test, + // but we replace these later on (apply_psf_model) + // set source->model based on best psf model fit + for (int i = 0; i < test->sources->n; i++) { + psSource *source = test->sources->data[i]; + // drop masked sources from PSFSTAR list + if (test->mask->data.U8[i]) { + source->type = PS_SOURCE_OTHER; + source->modelPSF = NULL; + } else { + // source->modelPSF = psMemCopy(test->modelPSF->data[i]); + source->modelPSF = test->modelPSF->data[i]; + } + } + # endif + + pmPSF *psf = psMemCopy(test->psf); + psFree (tests); + return (psf); +} Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 4946) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 4946) @@ -0,0 +1,114 @@ +# include "psphot.h" + +bool psphotFitGalaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) +{ + bool status, goodfit; + float x; + float y; + float sky; + int Nfit = 0; + int Nfail = 0; + int Niter = 0; + + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + + float GAL_MIN_SN = psMetadataLookupF32 (&status, config, "GAL_MIN_SN"); + float GAL_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA"); + float GAL_FIT_PADDING = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING"); + float GAL_MOMENTS_RAD = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS"); + + float FLUX_LIMIT = GAL_FIT_NSIGMA * skyStats->sampleStdev; + + char *modelName = psMetadataLookupPtr (&status, config, "GAL_MODEL"); + psModelType modelType = psModelSetType (modelName); + psModelRadius modelRadius = psModelRadius_GetFunction (modelType); + + psTimerStart ("psphot"); + + for (int i = 0; i < sources->n; i++) { + psSource *source = sources->data[i]; + + // sources which should not be fitted + // skip all valid stars + if (source->type == PS_SOURCE_PSFSTAR) continue; + if (source->type == PS_SOURCE_SATSTAR) continue; + if (source->type == PS_SOURCE_GOODSTAR) continue; + + // skip all likely defects + if (source->type == PS_SOURCE_DEFECT) continue; + if (source->type == PS_SOURCE_SATURATED) continue; + + // skip poorly fitted stars + if (source->type == PS_SOURCE_FAINTSTAR) continue; + if (source->type == PS_SOURCE_POOR_FIT_PSF) continue; + + // XXX EAM when do we pick these up again? + if (source->moments->SN < GAL_MIN_SN) { + source->type = PS_SOURCE_FAINT_GALAXY; // better choice? + continue; + } + + // recalculate the source moments using the larger galaxy moments radius + status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD); + if (!status) { + source->type = PS_SOURCE_DROP_GALAXY; // better choice? + continue; + } + + // use the source moments, etc to guess basic model parameters + psModel *model = pmSourceModelGuess (source, modelType); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // set the fit radius based on the object flux limit and the model + // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median) + model->radius = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING; + if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius"); + + // check if we need to redefine the pixels + // XXX EAM : a better test would examine the source pixels + if (model->radius > OUTER_RADIUS) { + // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, imdata, x, y, model->radius); + } + + // fit FLT (not PSF) model (set/unset the pixel mask) + psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP); + status = pmSourceFitModel_EAM (source, model, false); + psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP); + if (!status) { + // if the fit fails, we need to change the classification + psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + source->type = PS_SOURCE_FAIL_FIT_GAL; // better choice? + source->modelFLT = model; + Nfail ++; + continue; + } + + // check if model fit is acceptable + goodfit = pmModelFitStatus (model); + if (!goodfit) { + // if the fit fails, we need to change the classification + psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius); + source->type = PS_SOURCE_POOR_FIT_GAL; // better choice? + source->modelFLT = model; + Nfail ++; + continue; + } + + // accept the model + source->type = PS_SOURCE_GALAXY; + source->modelFLT = model; + Niter += model[0].nIter; + Nfit++; + + // subtract object, leave local sky + sky = model->params->data.F32[0]; + model->params->data.F32[0] = 0; + pmSourceSubModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + } + psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter); + return (true); +} Index: /trunk/psphot/src/psphotImageStats.c =================================================================== --- /trunk/psphot/src/psphotImageStats.c (revision 4946) +++ /trunk/psphot/src/psphotImageStats.c (revision 4946) @@ -0,0 +1,56 @@ +# include "psphot.h" + +psStats *psphotImageStats (psImageData *imdata, psMetadata *config) +{ + psStats *stats = NULL; + psStats *sky = NULL; + + psImage *image = imdata->image; + + // get image stats on a subset of the image (random 1e5 pts) + psTimerStart ("psphot"); + + // this section fixes a bug in the psImageStats definition: + // it should limit the number of elements used for the initial + // pass on the stats + { + bool status = false; + int Npixels = image->numRows*image->numCols; + int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "IMSTATS_NPIX")); + psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32); + psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); + + // choose Nsubset points between 0 and Nx*Ny, convert to coords + Npixels = image->numRows*image->numCols; + for (int i = 0; i < Nsubset; i++) { + double frnd = psRandomUniform (rnd); + int pixel = Npixels * frnd; + int ix = pixel / image->numCols; + int iy = pixel % image->numCols; + subset->data.F32[i] = image->data.F32[iy][ix]; + } + + // this should use ROBUST not SAMPLE median + // robust median is broken in pslib (0.5) + stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + stats = psVectorStats (stats, subset, NULL, NULL, 0); // use a mask? + psFree (subset); + psFree (rnd); + } + + // we store these values in mean,stdev + bool status = false; + float NOISE = psMetadataLookupF32 (&status, config, "RDNOISE"); + float GAIN = psMetadataLookupF32 (&status, config, "GAIN"); + + // convert instrumental background to poisson stats + sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + sky->sampleMean = stats->sampleMedian; + sky->sampleStdev = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE)); + + psLogMsg ("psphot", 3, "stats: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev); + + psFree (stats); + return (sky); +} Index: /trunk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 4946) +++ /trunk/psphot/src/psphotMarkPSF.c (revision 4946) @@ -0,0 +1,96 @@ +# include "psphot.h" + +// identify objects consistent with PSF shape/magnitude distribution +// we expect dparams[4],dparams[5] to have a scatter of: +// sigma_x / (S/N) * sqrt(2) +// 1 / (params[4],params[5])*(S/N) + +// sigma_x : 1 / SX +// dsx : 1 / (SX * SN) +// dsx_o = hypot (1/(SX*SN), MIN_DSX) + +// any objects which is consistent with the PSF should have +// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag) +// this includes a minimum buffer (DS) for the brighter objects + +// saturated stars should fall outside (but are already IDed) +// galaxies should be larger, cosmic rays smaller, but need to test? +// we also reject objects with S/N too low or ChiSquare to high + +// any object which meets the criterion is marked as PS_SOURCE_BRIGHTSTAR + +// floor for DS value +// XXX EAM : add to configuration? +# define MIN_DS 0.01 + +bool psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE) +{ + int keep; + float dSX, dSY, SX, SY, SN; + float nSx, nSy, Chi; + + if (source->modelPSF == NULL) return (false); + + // if object has fitted peak above saturation, label as SATSTAR + // remember: fit does not use saturated pixels (masked) + if (source->modelPSF->params->data.F32[1] >= SATURATE) { + if (source->type == PS_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n"); + } + source->type = PS_SOURCE_SATSTAR; + return (true); + } + if (source->type == PS_SOURCE_SATSTAR) { + psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); + source->type = PS_SOURCE_GOODSTAR; + } + + SN = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]; + SX = source->modelPSF->params->data.F32[4]; + SY = source->modelPSF->params->data.F32[5]; + dSX = source->modelPSF->dparams->data.F32[4]; + dSY = source->modelPSF->dparams->data.F32[5]; + Chi = source->modelPSF->chisq / source->modelPSF->nDOF; + + // swing of sigma_x,y in sigmas + nSx = dSX / hypot (MIN_DS, 1 / (SX * SN)); + nSy = dSY / hypot (MIN_DS, 1 / (SY * SN)); + + // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[] + keep = TRUE; + keep &= (fabs(nSx) < shapeNsigma); + keep &= (fabs(nSy) < shapeNsigma); + keep &= (SN > minSN); + keep &= (Chi < maxChi); + if (keep) { + if (source->type == PS_SOURCE_PSFSTAR) return (true); + source->type = PS_SOURCE_GOODSTAR; + return (true); + } + + if (source->type == PS_SOURCE_PSFSTAR) { + psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n"); + } + + // object appears to be small, suspected defect + if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) { + source->type = PS_SOURCE_DEFECT; + return (false); + } + + // object appears to be large, suspected galaxy + if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) { + source->type = PS_SOURCE_GALAXY; + return (false); + } + + // object appears to be extremely faint: what is this? + if (SN < minSN) { + source->type = PS_SOURCE_FAINTSTAR; + return (false); + } + + // these are pooly fitted, probable stars near other stars? + source->type = PS_SOURCE_POOR_FIT_PSF; + return (false); +} Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 4946) +++ /trunk/psphot/src/psphotOutput.c (revision 4946) @@ -0,0 +1,643 @@ +# include "psphot.h" + +// output functions: we have several fixed modes (sx, obj, cmp) +void psphotOutput (psImageData *imdata, psMetadata *config, psArray *sources) { + + bool status; + + char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE"); + char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE"); + char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE"); + + if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage); + + if (outputFile == NULL) { + psLogMsg ("output", 3, "no data output file selected"); + return; + } + if (outputMode == NULL) { + psLogMsg ("output", 3, "no data output mode selected"); + return; + } + if (!strcasecmp (outputMode, "TEXT")) { + pmSourcesWriteText (imdata, outputFile, sources); + return; + } + + if (!strcasecmp (outputMode, "OBJ")) { + pmSourcesWriteOBJ (imdata, outputFile, sources); + return; + } + + if (!strcasecmp (outputMode, "SX")) { + pmSourcesWriteSX (imdata, outputFile, sources); + return; + } + + if (!strcasecmp (outputMode, "CMP")) { + pmSourcesWriteCMP (imdata, outputFile, sources); + return; + } + + if (!strcasecmp (outputMode, "CMF")) { + pmSourcesWriteCMF (imdata, outputFile, sources); + return; + } + psAbort ("psphot", "unknown output mode %s", outputMode); +} + +bool pmSourcesWriteText (psImageData *imdata, char *filename, psArray *sources) { + + char *name = (char *) psAlloc (strlen(filename) + 10); + + sprintf (name, "%s.psf.dat", filename); + psModelWritePSFs (sources, name); + + sprintf (name, "%s.flt.dat", filename); + psModelWriteFLTs (sources, name); + + sprintf (name, "%s.nul.dat", filename); + psModelWriteNULLs (sources, name); + + sprintf (name, "%s.mnt.dat", filename); + psMomentsWriteText (sources, name); + + psFree (name); + return true; +} + +// dophot-style output list with fixed line width +bool pmSourcesWriteOBJ (psImageData *imdata, char *filename, psArray *sources) { + + int i, type, status; + psModel *model; + psF32 *PAR, *dPAR; + float sky, dmag, apMag, fitMag; + + 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 first + for (i = 0; i < sources->n; i++) { + psSource *source = (psSource *) sources->data[i]; + + // use the correct model + switch (source->type) { + + // use PSF model with stars + case PS_SOURCE_PSFSTAR: + case PS_SOURCE_SATSTAR: + case PS_SOURCE_GOODSTAR: + model = source->modelPSF; + break; + + // use FLT model with galaxies + case PS_SOURCE_GALAXY: + model = source->modelFLT; + break; + + // skip defects and poorly fitted stars or galaxies + case PS_SOURCE_DEFECT: + case PS_SOURCE_SATURATED: + case PS_SOURCE_FAINTSTAR: + case PS_SOURCE_POOR_FIT_PSF: + case PS_SOURCE_POOR_FIT_GAL: + case PS_SOURCE_FAIL_FIT_GAL: + default: + model = NULL; + } + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // save local sky for later + sky = model->params->data.F32[0]; + + // replace source flux + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceAddModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // measure object photometry + status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // subtract object, leave local sky + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceSubModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + if (status == FALSE) continue; + + dmag = dPAR[1] / PAR[1]; + type = pmSourceDophotType (source); + psLineInit (line); + psLineAdd (line, "%3d", source->type); + psLineAdd (line, "%8.2f", PAR[2]); + psLineAdd (line, "%8.2f", PAR[3]); + psLineAdd (line, "%8.3f", 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", 32.0); + psLineAdd (line, "%8.3f", apMag); + psLineAdd (line, "%8.2f\n", apMag - fitMag); + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// elixir/sextractor-style output list with fixed line width +bool pmSourcesWriteSX (psImageData *imdata, char *filename, psArray *sources) { + + int i, status, flags; + psModel *model; + psF32 *PAR, *dPAR; + float sky, dmag, apMag, fitMag; + + psLine *line = psLineAlloc (110); // 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 first + for (i = 0; i < sources->n; i++) { + psSource *source = (psSource *) sources->data[i]; + flags = 0; + + // use the correct model + switch (source->type) { + + // use PSF model with stars + case PS_SOURCE_SATSTAR: + flags = 4; + model = source->modelPSF; + break; + + case PS_SOURCE_PSFSTAR: + case PS_SOURCE_GOODSTAR: + model = source->modelPSF; + break; + + // use FLT model with galaxies + case PS_SOURCE_GALAXY: + model = source->modelFLT; + break; + + // skip defects and poorly fitted stars or galaxies + case PS_SOURCE_DEFECT: + case PS_SOURCE_SATURATED: + case PS_SOURCE_FAINTSTAR: + case PS_SOURCE_POOR_FIT_PSF: + case PS_SOURCE_POOR_FIT_GAL: + case PS_SOURCE_FAIL_FIT_GAL: + default: + model = NULL; + } + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // save local sky for later + sky = model->params->data.F32[0]; + + // replace source flux + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceAddModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // measure object photometry + status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // subtract object, leave local sky + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceSubModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + if (status == FALSE) continue; + + // XXX EAM : fix this to match sextractor types? + // type = pmSourceSextractType (source); + + 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", 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", 32.0); // should be MAG_ISO + psLineAdd (line, "%9.4f", 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 (psImageData *imdata, char *filename, psArray *sources) { + + int i, type, status; + psModel *model; + psF32 *PAR, *dPAR; + float sky, dmag, apMag, fitMag, lsky, Theta; + + // create file, write-out header + unlink (filename); + psFits *fits = psFitsAlloc (filename); + psFitsWriteHeader (imdata->header, fits); + psFree (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 (104); // 104 is dophot-defined line length + + // write sources with models first + for (i = 0; i < sources->n; i++) { + psSource *source = (psSource *) sources->data[i]; + + // use the correct model + switch (source->type) { + + // use PSF model with stars + case PS_SOURCE_PSFSTAR: + case PS_SOURCE_SATSTAR: + case PS_SOURCE_GOODSTAR: + model = source->modelPSF; + break; + + // use FLT model with galaxies + case PS_SOURCE_GALAXY: + model = source->modelFLT; + break; + + // skip defects and poorly fitted stars or galaxies + case PS_SOURCE_DEFECT: + case PS_SOURCE_SATURATED: + case PS_SOURCE_FAINTSTAR: + case PS_SOURCE_POOR_FIT_PSF: + case PS_SOURCE_POOR_FIT_GAL: + case PS_SOURCE_FAIL_FIT_GAL: + default: + model = NULL; + } + if (model == NULL) continue; + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + + // save local sky for later + sky = model->params->data.F32[0]; + + // XXX EAM : get Theta from PAR[4,5,6] + Theta = 0.0; + + // replace source flux + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceAddModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // measure object photometry + status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask); + + // subtract object, leave local sky + // XXX EAM : need to add 'sky' boolean option here + model->params->data.F32[0] = 0; + pmSourceSubModel (source->pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + if (status == FALSE) continue; + + 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 ", fitMag + 25.0); + psLineAdd (line, "%03d ", (int)(1000*dmag)); + psLineAdd (line, "%2d ", type); + psLineAdd (line, "%3.1f ", lsky); + psLineAdd (line, "%6.3f ", 32.0); // should be 'Mgal + psLineAdd (line, "%6.3f ", apMag + 25.0); + 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", Theta); // should be theta + fwrite (line->line, 1, line->Nline, f); + } + fclose (f); + return true; +} + +// elixir-style FITS table output (header + table in 1st extension) +bool pmSourcesWriteCMF (psImageData *imdata, char *filename, psArray *sources) { + + // write imdata->header to file + + // write FITS table to file (use autocode tools) + + return true; +} + +/***** Text Output Methods *****/ + +// write the peaks to an output file +bool psPeaksWriteText (psArray *sources, char *filename) { + + int i; + FILE *f; + psPeak *peak; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("psPeaksWriteText", 3, "can't open output file for peaks%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + peak = (psPeak *) sources->data[i]; + if (peak == NULL) continue; + fprintf (f, "%5d %5d %7.1f\n", + peak->x, peak->y, peak->counts); + } + fclose (f); + return true; +} + +// write the moments to an output file +bool psMomentsWriteText (psArray *sources, char *filename) { + + int i; + FILE *f; + psSource *source; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("psMomentsWriteText", 3, "can't open output file for moments%s\n", filename); + return false; + } + + for (i = 0; i < sources->n; i++) { + source = (psSource *) sources->data[i]; + if (source->moments == NULL) continue; + 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); + return true; +} + +// dump the sources to an output file +bool psModelWritePSFs (psArray *sources, char *filename) { + + double dP, flux; + int i, j; + FILE *f; + psVector *params; + psVector *dparams; + psModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("psModelWritePSFs", 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++) { + psSource *source = (psSource *) sources->data[i]; + model = (psModel *) source->modelPSF; + if (model == NULL) continue; + + // valid source types for this function + if (source->type == PS_SOURCE_SATSTAR) goto valid; + if (source->type == PS_SOURCE_PSFSTAR) goto valid; + if (source->type == PS_SOURCE_GOODSTAR) goto valid; + continue; + + valid: + params = model->params; + dparams = model->dparams; + + dP = 0; + dP += PS_SQR(dparams[0].data.F32[2]); + dP += PS_SQR(dparams[0].data.F32[3]); + dP = sqrt (dP); + + psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type); + flux = modelFluxFunc (params); + + fprintf (f, "%7.1f %7.1f %5.1f %7.3f %7.4f %7.4f ", + params[0].data.F32[2], + params[0].data.F32[3], + params[0].data.F32[0], + -2.5*log10(flux), + (dparams[0].data.F32[1]/params[0].data.F32[1]), + dP); + for (j = 0; j < model->params->n - 4; j++) { + fprintf (f, "%9.6f ", params[0].data.F32[j+4]); + } + fprintf (f, " : "); + for (j = 0; j < model->params->n - 4; j++) { + fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]); + } + fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", + source[0].type, + 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 psModelWriteFLTs (psArray *sources, char *filename) { + + double dP; + int i, j; + FILE *f; + psVector *params; + psVector *dparams; + psModel *model; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("psModelWriteFLTs", 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++) { + psSource *source = (psSource *) sources->data[i]; + model = (psModel *) source->modelFLT; + if (model == NULL) continue; + if (source->type == PS_SOURCE_GALAXY) goto valid; + continue; + + valid: + params = model->params; + dparams = model->dparams; + + // XXX these are hardwired for SGAUSS : this should be pushed into the + // model functions as an abstract function + dP = 0; + dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]); + dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]); + dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]); + dP = sqrt (dP); + + fprintf (f, "%7.1f %7.1f %5.1f %7.3f %7.4f %7.4f ", + params[0].data.F32[2], + params[0].data.F32[3], + params[0].data.F32[0], + -2.5*log10(params[0].data.F32[1]), + (dparams[0].data.F32[1]/params[0].data.F32[1]), + dP); + + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", params[0].data.F32[j]); + } + fprintf (f, " : "); + for (j = 4; j < model->params->n; j++) { + fprintf (f, "%9.6f ", dparams[0].data.F32[j]); + } + fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", + source[0].type, 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 psModelWriteNULLs (psArray *sources, char *filename) +{ + + int i; + FILE *f; + + f = fopen (filename, "w"); + if (f == NULL) { + psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename); + return false; + } + + psMoments *empty = pmMomentsAlloc (); + + // write sources with models first + for (i = 0; i < sources->n; i++) { + psSource *source = (psSource *) sources->data[i]; + + // skip these sources (in PSF or FLT) + if (source->type == PS_SOURCE_GALAXY) continue; + if (source->type == PS_SOURCE_PSFSTAR) continue; + if (source->type == PS_SOURCE_SATSTAR) continue; + if (source->type == PS_SOURCE_GOODSTAR) continue; + + if (source->moments == NULL) { + 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, + empty->x, empty->y, + empty->Sx, empty->Sy, + empty->Sum, empty->Peak, + empty->Sky, empty->SN, + empty->nPixels, source->type); + } else { + 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); + return true; +} + +// translations between psphot object types and dophot object types +int pmSourceDophotType (psSource *source) { + + switch (source->type) { + + case PS_SOURCE_DEFECT: + case PS_SOURCE_SATURATED: + return (8); + + case PS_SOURCE_SATSTAR: + return (10); + + case PS_SOURCE_PSFSTAR: + case PS_SOURCE_GOODSTAR: + return (1); + + case PS_SOURCE_POOR_FIT_PSF: + return (7); + + case PS_SOURCE_FAIL_FIT_PSF: + case PS_SOURCE_FAINTSTAR: + return (4); + + case PS_SOURCE_GALAXY: + case PS_SOURCE_FAINT_GALAXY: + case PS_SOURCE_DROP_GALAXY: + case PS_SOURCE_FAIL_FIT_GAL: + case PS_SOURCE_POOR_FIT_GAL: + return (2); + + case PS_SOURCE_OTHER: + default: + return (0); + } + return (0); +} + +int psphotSaveImage (psMetadata *header, psImage *image, char *filename) { + + unlink (filename); + psFits *fits = psFitsAlloc (filename); + psFitsWriteImage (fits, header, image, 0, NULL); + psFree (fits); + return (TRUE); +} Index: /trunk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 4946) +++ /trunk/psphot/src/psphotSetup.c (revision 4946) @@ -0,0 +1,96 @@ +# include "psphot.h" + +// load the image +// load or construct the mask +// apply X/Y border info +// apply the SATURATE mask? +// load or construct the noise + +psImageData *psphotSetup (psMetadata *config) +{ + psMetadata *header = NULL; + psImage *image = NULL; + psImage *noise = NULL; + psImage *mask = NULL; + psRegion region = {0,0,0,0}; // a region representing the entire array + bool status; + + psTimerStart ("psphot"); + + // setup header template, specify COMMENT and HISTORY as MULTI + // XXX EAM - shouldn't this be part of the basic FITS concept? + header = psMetadataAlloc (); + psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL); + psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL); + + // read header and image data from INPUT + char *input = psMetadataLookupPtr (&status, config, "IMAGE"); + psFits *file = psFitsAlloc (input); + header = psFitsReadHeader (header, file); + image = psFitsReadImage (NULL, file, region, 0); + psFree (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, header, "RDNOISE"); + float GAIN = pmConfigLookupF32 (&status, config, header, "GAIN"); + psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); + psMetadataAdd (config, PS_LIST_TAIL, "GAIN", PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN); + + // load the noise image if it is supplied, otherwise build from input + char *noiseName = psMetadataLookupPtr (&status, config, "NOISE"); + if (status == true) { + file = psFitsAlloc (noiseName); + noise = psFitsReadImage (NULL, file, region, 0); + psFree (file); + // psFree (noiseName); XXX - see psFree (input) + } else { + // build the noise image from the input image, RDNOISE, and GAIN + noise = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32); + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + noise->data.F32[iy][ix] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN); + } + } + } + + // load the mask if specified, otherwise construct + char *maskName = psMetadataLookupPtr (&status, config, "MASK"); + if (status == true) { + // XXX EAM require / convert mask to psU8? + file = psFitsAlloc (maskName); + mask = psFitsReadImage (NULL, file, region, 0); + psFree (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 = psRegionAlloc (XMIN, XMAX, YMIN, YMAX); + keep = psRegionForImage (keep, image, keep); + psImageKeepRegion (mask, keep, OR, PSPHOT_MASK_INVALID); + psFree (keep); + + // mask the saturated pixels + // XXX EAM does the mask need to grow? + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + for (int i = 0; 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; + } + } + } + + psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); + + // save the image data & return it + psImageData *imdata = psImageDataAlloc(image, noise, mask, header); + return (imdata); +} Index: /trunk/psphot/src/psphotSortBySN.c =================================================================== --- /trunk/psphot/src/psphotSortBySN.c (revision 4946) +++ /trunk/psphot/src/psphotSortBySN.c (revision 4946) @@ -0,0 +1,19 @@ +# include "psphot.h" + +// sort by SN (descending) +int psphotSortBySN (const void **a, const void **b) +{ + psSource *A = *(psSource **)a; + psSource *B = *(psSource **)b; + + psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN; + psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN; + if (isnan (fA)) fA = 0; + if (isnan (fB)) fB = 0; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (-1); + if (diff < FLT_EPSILON) return (+1); + return (0); +} + Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 4946) +++ /trunk/psphot/src/psphotSourceStats.c (revision 4946) @@ -0,0 +1,54 @@ +# include "psphot.h" + +psArray *psphotSourceStats (psImageData *imdata, psMetadata *config, psArray *peaks) +{ + bool status = false; + psArray *sources = NULL; + + psTimerStart ("psphot"); + + // determine properties (sky, moments) of initial sources + float INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS"); + float OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS"); + + sources = psArrayAlloc (peaks->n); + sources->n = 0; + + for (int i = 0; i < peaks->n; i++) { + + // create a new source, add peak + psSource *source = pmSourceAlloc(); + source->peak = (psPeak *)psMemCopy(peaks->data[i]); + + // allocate image, noise, mask arrays for each peak (square of radius OUTER) + // XXX EAM : this has 12 leaks (v.5) + pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER); + + // measure a local sky value + // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER); + if (!status) { + psFree (source); + continue; + } + + // measure basic source moments + status = pmSourceMoments (source, RADIUS); + if (!status) { + psFree (source); + continue; + } + + // add to the source array + psArrayAdd (sources, 100, source); + psFree (source); + } + + psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", 3, "%d sources\n", sources->n); + + return (sources); +} + +// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky) Index: /trunk/psphot/src/psphotSubtractPSF.c =================================================================== --- /trunk/psphot/src/psphotSubtractPSF.c (revision 4946) +++ /trunk/psphot/src/psphotSubtractPSF.c (revision 4946) @@ -0,0 +1,36 @@ +# include "psphot.h" + +// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) +// need to control application of this with some thresholding + +bool psphotSubtractPSF (psSource *source) +{ + float sky; + psModel *model; + psImage *pixels; + + // only subtract successful fits + if (source->type == PS_SOURCE_SATSTAR) goto valid; + if (source->type == PS_SOURCE_PSFSTAR) goto valid; + if (source->type == PS_SOURCE_GOODSTAR) goto valid; + return (false); + +valid: + model = source->modelPSF; + if (model == NULL) { + psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y); + return (false); + } + + pixels = source->pixels; + + // subtract object, leave local sky + sky = model->params->data.F32[0]; + model->params->data.F32[0] = 0; + pmSourceSubModel (pixels, source->mask, model, false); + model->params->data.F32[0] = sky; + + // XXX EAM : amplify the noise matrix a la dophot? + + return (true); +} Index: unk/psphot/src/setup.c =================================================================== --- /trunk/psphot/src/setup.c (revision 4945) +++ (revision ) @@ -1,108 +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 noise - -psImageData *setup (psMetadata *config) -{ - psMetadata *header = NULL; - psImage *image = NULL; - psImage *noise = NULL; - psImage *mask = NULL; - psRegion region = {0,0,0,0}; // a region representing the entire array - bool status; - - psTimerStart ("psphot"); - - // setup header template, specify COMMENT and HISTORY as MULTI - header = psMetadataAlloc (); - psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL); - psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL); - - // read header and image data from INPUT - char *input = psMetadataLookupPtr (&status, config, "INPUT"); - psFits *file = psFitsAlloc (input); - header = psFitsReadHeader (header, file); - image = psFitsReadImage (NULL, file, region, 0); - psFree (file); - // psFree (input); - // XXX - 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, header, "RDNOISE"); - float GAIN = pmConfigLookupF32 (&status, config, header, "GAIN"); - psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE); - psMetadataAdd (config, PS_LIST_TAIL, "GAIN", PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN); - - // load the noise image if it is supplied, otherwise build from input - char *noiseName = psMetadataLookupPtr (&status, config, "NOISE"); - if (status == true) { - file = psFitsAlloc (noiseName); - noise = psFitsReadImage (NULL, file, region, 0); - psFree (file); - // psFree (noiseName); XXX - see psFree (input) - } else { - - // XXX EAM it looks like psBinaryOp is broken for: image op scalar - // this loop performs the operation by hand - noise = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32); - for (int iy = 0; iy < image->numRows; iy++) { - for (int ix = 0; ix < image->numCols; ix++) { - noise->data.F32[iy][ix] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN); - } - } - - // this uses psBinaryOp - # if (0) - psScalar *value; - - value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64); - noise = (psImage *) psBinaryOp (NULL, image, "/", value); - psFree (value); - - value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64); - noise = (psImage *) psBinaryOp (noise, noise, "+", value); - psFree (value); - # endif - } - - // load the mask if specified, otherwise construct - char *maskName = psMetadataLookupPtr (&status, config, "MASK"); - if (status == true) { - file = psFitsAlloc (maskName); - mask = psFitsReadImage (NULL, file, region, 0); - psFree (file); - // psFree (maskName); XXX - see psFree (input) - // XXX require / convert mask to psU8? - } else { - mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8); - } - - 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 = psRegionAlloc (XMIN, XMAX, YMIN, YMAX); - keep = psRegionForImage (keep, image, keep); - psImageKeepRegion (mask, keep, OR, 0x01); - psFree (keep); - - // mask the saturated pixels - float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE"); - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - if (image->data.F32[i][j] >= SATURATE) { - mask->data.U8[i][j] |= 0x02; - } - } - } - - psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot")); - - // save the image data & return it - psImageData *imdata = psImageDataAlloc(image, noise, mask, header); - return (imdata); -} Index: unk/psphot/src/source_moments.c =================================================================== --- /trunk/psphot/src/source_moments.c (revision 4945) +++ (revision ) @@ -1,53 +1,0 @@ -# include "psphot.h" - -psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *peaks) -{ - bool status = false; - psArray *sources = NULL; - - psMemoryId id = psMemGetId (); - - psTimerStart ("psphot"); - - // determine properties (sky, moments) of initial sources - float INNER = psMetadataLookupF32 (&status, config, "INNER_RADIUS"); - float OUTER = psMetadataLookupF32 (&status, config, "OUTER_RADIUS"); - float RADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS"); - - sources = psArrayAlloc (peaks->n); - sources->n = 0; - - for (int i = 0; i < peaks->n; i++) { - - // create a new source, add peak - psSource *source = pmSourceAlloc(); - source->peak = (psPeak *)psMemCopy(peaks->data[i]); - - // allocate image, noise, mask arrays for each peak (square of radius OUTER) - // XXX - this has 12 leaks (v.5) - pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER); - - // this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER); - if (!status) { - psFree (source); - continue; - } - - status = pmSourceMoments (source, RADIUS); - if (!status) { - psFree (source); - continue; - } - - psArrayAdd (sources, 100, source); - psFree (source); - } - - psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot")); - psLogMsg ("psphot", 3, "%d sources\n", sources->n); - - return (sources); -} - -// filter out bad peaks (eg, no valid pixels available for sky) Index: unk/psphot/src/subtract_psf_source.c =================================================================== --- /trunk/psphot/src/subtract_psf_source.c (revision 4945) +++ (revision ) @@ -1,36 +1,0 @@ -# include "psphot.h" - -// subtract model fits from image (only PSFSTAR && BRIGHTSTAR) -// need to control application of this with some thresholding - -bool subtract_psf_source (psSource *source) -{ - float sky; - psModel *model; - psImage *pixels; - - // only subtract successful fits - if (source->type == PS_SOURCE_SATSTAR) goto valid; - if (source->type == PS_SOURCE_PSFSTAR) goto valid; - if (source->type == PS_SOURCE_GOODSTAR) goto valid; - return (false); - -valid: - model = source->modelPSF; - if (model == NULL) { - psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y); - return (false); - } - - pixels = source->pixels; - - // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - - // XXX adjust (enhance) the noise matrix a la dophot? - - return (true); -} Index: unk/psphot/src/test_psf_scatter.c =================================================================== --- /trunk/psphot/src/test_psf_scatter.c (revision 4945) +++ (revision ) @@ -1,82 +1,0 @@ -# include "psphot.h" - -// measure the distribution of dparams[4,5] as a function of magnitude -// use this as a test comparison with the theoretical distribution - -bool test_psf_scatter (psArray *sources) -{ - // analysis of PSF shape distribution - // select the PSF sources - int Npt = 0; - psVector *deltaSX = psVectorAlloc (100, PS_TYPE_F64); - psVector *deltaSY = psVectorAlloc (100, PS_TYPE_F64); - psVector *logSN = psVectorAlloc (100, PS_TYPE_F64); - for (int i = 0; i < sources->n; i++) { - psSource *source = sources->data[i]; - if (source->type != PS_SOURCE_PSFSTAR) continue; - if (source->modelPSF == NULL) continue; - deltaSX->data.F64[Npt] = source->modelPSF->dparams->data.F32[4]; - deltaSY->data.F64[Npt] = source->modelPSF->dparams->data.F32[5]; - logSN ->data.F64[Npt] = log10(source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]); - Npt ++; - if (Npt == deltaSX->nalloc) { - deltaSX = psVectorRealloc (deltaSX, deltaSX->nalloc + 100); - deltaSY = psVectorRealloc (deltaSY, deltaSY->nalloc + 100); - logSN = psVectorRealloc (logSN, logSN ->nalloc + 100); - } - } - if (Npt < 45) psAbort ("PSF analysis", "too few PSF stars for scatter analysis"); - logSN->n = deltaSX->n = deltaSY->n = Npt; - - // use a min number of bins or a min number of stars per bin? - // determine dynamic range of magnitudes - psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); - psVectorStats (stats, logSN, NULL, NULL, 0); - float snMin = stats->min; - float snMax = stats->max; - float snDel = 25*(snMax - snMin) / logSN->n; - int nBins = 1 + (int)((snMax - snMin) / snDel); - - // generate a log-Sx, mag line and fit for appropriate mag bins - stats->options = PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV; - psVector *snBin = psVectorAlloc (nBins, PS_TYPE_F64); - psVector *dsxBin = psVectorAlloc (nBins, PS_TYPE_F64); - psVector *dsyBin = psVectorAlloc (nBins, PS_TYPE_F64); - psVector *tsx = psVectorAlloc (logSN->n, PS_TYPE_F64); - psVector *tsy = psVectorAlloc (logSN->n, PS_TYPE_F64); - int bin = 0; - for (float sn = snMin; sn <= snMax - snDel; sn += snDel) { - tsx->n = tsy->n = 0; - for (int i = 0; i < logSN->n; i++) { - if (logSN->data.F64[i] < sn) continue; - if (logSN->data.F64[i] > sn + snDel) continue; - tsx->data.F64[tsx->n] = deltaSX->data.F64[i]; - tsy->data.F64[tsy->n] = deltaSY->data.F64[i]; - tsx->n++; - tsy->n++; - } - if (tsx->n < 5) continue; - snBin ->data.F64[bin] = sn + 0.5*snDel; - psVectorStats (stats, tsx, NULL, NULL, 0); - dsxBin->data.F64[bin] = log10(stats->clippedStdev); - psVectorStats (stats, tsy, NULL, NULL, 0); - dsyBin->data.F64[bin] = log10(stats->clippedStdev); - fprintf (stderr, "bin: %d, sn: %f, ndsx: %d, ndsy: %d, dsx: %f, dsy: %f\n", - bin, sn + 0.5*snDel, tsx->n, tsy->n, dsxBin->data.F64[bin], dsyBin->data.F64[bin]); - bin ++; - } - dsxBin->n = bin; - dsyBin->n = bin; - snBin ->n = bin; - - // XXX these APIs should get changed - psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD); - psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL); - psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD); - psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL); - // these should have a slope of 1.0 - psPolynomial1DDump (dsxLine); - psPolynomial1DDump (dsyLine); - - return (true); -}