Changeset 4582

Timestamp:

Jul 20, 2005, 4:19:11 AM (21 years ago)

Author:

eugene

Message:

working on the stability of the galaxy models

Location:

trunk/psphot/src

Files:

• LocalSky.c (modified) (1 diff)
• apply_psf_model.c (modified) (1 diff)
• basic_classes.c (modified) (1 diff)
• fit_galaxies.c (modified) (3 diffs)
• onesource.c (modified) (2 diffs)
• psphot.h (modified) (2 diffs)

trunk/psphot/src/LocalSky.c

@@ -58 +58 @@
 }
 
+# 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,
+                          float SOFT)
+{
+    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;
+
+
+    // XXX EAM : is it necessary for the mask & noise to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    int nParams = PSF ? params->n - 4 : params->n;
+    psF32 Xo = params->data.F32[2];
+    psF32 Yo = params->data.F32[3];
+
+    // 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];
+
+                // XXX just for test purposes, use the raw radius.  a better
+                // choice is to ask what is the value of z and scale by that
+                Ro = hypot ((Xo - coord->data.F32[0]), (Yo - coord->data.F32[1]));
+
+                // XXX enhance the noise (doubled if R = 10 pix, etc. note the square)
+                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * (1 + PS_SQR(Ro/SOFT)));
+                // 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, params->n, PS_TYPE_F64);
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2(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);
+}
+

trunk/psphot/src/apply_psf_model.c

@@ -21 +21 @@
     float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
     float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
-    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
     float OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
 

trunk/psphot/src/basic_classes.c

@@ -11 +11 @@
     // make this optional
     // DumpMoments (sources, "moments.dat");
+    return (true);
 }

trunk/psphot/src/fit_galaxies.c

@@ -14 +14 @@
 
     float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
-    float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+    // float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
     float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
     float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
@@ -22 +22 @@
     float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
 
-    psModelType   modelType   = psModelSetType ("PS_MODEL_RGAUSS");
+    psModelType   modelType   = psModelSetType ("PS_MODEL_SGAUSS");
     psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
 
@@ -34 +34 @@
 
         // recalculate the source moments using the galaxy radius (larger)
-        status = pmSourceMoments (source, MOMENT_R);
+        status = pmSourceMoments_EAM (source, MOMENT_R);
+        if (!status) {
+          fprintf (stderr, "invalid moments, skipping\n");
+          continue;
+        }
 
         // use the source moments, etc to guess basic model parameters

trunk/psphot/src/onesource.c

@@ -14 +14 @@
     float MRAD   = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
     float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+    float SOFT   = psMetadataLookupF32 (&status, config, "SOFT_RADIUS");
 
     psSource *source = pmSourceAlloc();
@@ -26 +27 @@
 
     psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
-    status = pmSourceFitModel (source, model, false);
+    status = pmSourceFitModel_EAM (source, model, false, SOFT);
 
     pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);

trunk/psphot/src/psphot.h

@@ -67 +67 @@
 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);
 
 // psphot utilities
@@ -77 +80 @@
 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, float SOFT);
+bool pmSourceMoments_EAM(psSource *source, psF32 radius);
 
 // fitsource utilities