Changeset 5993 for trunk/psphot/src/psphotApResid.c

Timestamp:

Jan 15, 2006, 8:30:16 AM (21 years ago)

Author:

eugene

Message:

API cleanup, removed old test files, some re-organization, changed FLT/GAL to EXT

File:

• trunk/psphot/src/psphotApResid.c (modified) (3 diffs)

trunk/psphot/src/psphotApResid.c

@@ -1 +1 @@
 # include "psphot.h"
-
-psPolynomial4D *psVectorChiClipFitPolynomial4D(
-    psPolynomial4D *poly,
-    psStats *stats,
-    const psVector *mask,
-    psMaskType maskValue,
-    const psVector *f,
-    const psVector *fErr,
-    const psVector *x,
-    const psVector *y,
-    const psVector *z,
-    const psVector *t)
-{
-    PS_ASSERT_POLY_NON_NULL(poly, NULL);
-    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
-    PS_ASSERT_PTR_NON_NULL(stats, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
-    if (mask != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
-        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
-    }
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(fErr, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
-
-    // clipping range defined by min and max and/or clipSigma
-    float minClipSigma;
-    float maxClipSigma;
-    if (isfinite(stats->max)) {
-        maxClipSigma = +fabs(stats->max);
-    } else {
-        maxClipSigma = +fabs(stats->clipSigma);
-    }
-    if (isfinite(stats->min)) {
-        minClipSigma = -fabs(stats->min);
-    } else {
-        minClipSigma = -fabs(stats->clipSigma);
-    }
-    psVector *fit   = NULL;
-    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
-
-    // eventual expansion: user supplies one of various stats option pairs,
-    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
-    // evaluate the clipping sigma
-    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
-    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-
-    for (int N = 0; N < stats->clipIter; N++) {
-        int Nkeep = 0;
-
-        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
-        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
-        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
-
-        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
-        psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
-
-        // set mask if pts are not valid
-        // we are masking out any point which is out of range
-        // recovery is not allowed with this scheme
-        for (int i = 0; i < resid->n; i++) {
-            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
-                continue;
-            }
-            float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev);
-            if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) {
-                if (mask != NULL) {
-                    mask->data.U8[i] |= 0x01;
-                }
-                continue;
-            }
-            if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) {
-                if (mask != NULL) {
-                    mask->data.U8[i] |= 0x01;
-                }
-                continue;
-            }
-            Nkeep ++;
-        }
-
-        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
-                 Nkeep, x->n);
-
-        stats->clippedNvalues = Nkeep;
-        psFree (fit);
-    }
-    // Free local temporary variables
-    psFree (resid);
-
-    if (poly == NULL) {
-        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
-        return(NULL);
-    }
-    return(poly);
-}
 
 // measure the aperture residual statistics
@@ -347 +239 @@
     psf->nApResid = residStats->clippedNvalues;
 
-    /*
-      (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
-      (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
-      (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
-    */
+    // save results for later output
+    psMetadataAdd (config, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
+    psMetadataAdd (config, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
+    psMetadataAdd (config, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
+    psMetadataAdd (config, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
+    psMetadataAdd (config, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->growth->apLoss);
+    psMetadataAdd (config, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
 
     psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
@@ -370 +264 @@
 }
 
+/*
+  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
+  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
+  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
+*/
+