Changeset 5048 for trunk/psphot/src/psModulesUtils.c

Timestamp:

Sep 13, 2005, 3:35:20 PM (21 years ago)

Author:

eugene

Message:

substantial work to match with psLib v.7

File:

• trunk/psphot/src/psModulesUtils.c (modified) (2 diffs)

trunk/psphot/src/psModulesUtils.c

@@ -2 +2 @@
 
 // extract config informatin from config data or from image header, if specified
+// XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme
 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
 
@@ -48 +49 @@
 }
 
-// XXX EAM : these are my alternate implementations of psModule functions
-
-// this sets and clears bit 0x80
-bool pmSourceLocalSky_EAM (pmSource *source,
-                           psStatsOptions statsOptions,
-                           psF32 Radius)
-{
-
-    psImage *image = source->pixels;
-    psImage *mask  = source->mask;
-    pmPeak *peak  = source->peak;
-    psRegion srcRegion;
-
-    // XXX EAM : psRegionXXX funcs need value not ptr
-
-    srcRegion = psRegionForSquare (peak->x, peak->y, Radius);
-    srcRegion = psRegionForImage (mask, &srcRegion);
-
-    psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED);
-    psStats *myStats = psStatsAlloc(statsOptions);
-    myStats = psImageStats(myStats, image, mask, 0xff);
-    psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED);
-
-    psF64 tmpF64;
-    p_psGetStatValue(myStats, &tmpF64);
-    psFree (myStats);
-
-    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 (pmSource *source,
-                             pmModel *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)
-
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-    pmModelLimits modelLimits = pmModelLimits_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, &params_min, &params_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
-pmSource data structure, along with the peak location, and determines the
-various moments associated with that peak.
- 
-Requires the following to have been created:
-    pmSource
-    pmSource->peak
-    pmSource->pixels
- 
-XXX: The peak calculations are done in image coords, not subImage coords.
- 
-XXX: mask values?
-
-XXX EAM : this version clips input pixels on S/N
-*****************************************************************************/
-# define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
-
-bool pmSourceMoments_EAM(pmSource *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 (pmModel *model) {