Changeset 4129

Timestamp:

Jun 6, 2005, 11:20:45 PM (21 years ago)

Author:

eugene

Message:

further basic concepts

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• .cvsignore (added)
• Makefile (modified) (1 diff)
• src/apply_psf_model.c (modified) (2 diffs)
• src/by_SN.c (modified) (2 diffs)
• src/choose_psf_model.c (modified) (1 diff)
• src/find_defects.c (added)
• src/mark_psf_source.c (modified) (1 diff)
• src/psPolynomials.c (modified) (2 diffs)
• src/psUtils.c (modified) (1 diff)
• src/psphot.c (modified) (4 diffs)
• src/psphot.h (modified) (4 diffs)
• src/pspsf.c (modified) (4 diffs)
• src/source_moments.c (modified) (1 diff)
• src/subtract_galaxies.c (modified) (1 diff)
• src/subtract_psf_source.c (modified) (2 diffs)
• src/test_psf_scatter.c (modified) (1 diff)

trunk/psphot/Makefile

@@ -32 +32 @@
 $(SRC)/mark_psf_sources.$(ARCH).o \
 $(SRC)/subtract_psf_sources.$(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)/subtract_galaxies.$(ARCH).o
+$(SRC)/subtract_galaxies.$(ARCH).o \
+$(SRC)/find_defects.$(ARCH).o
 
 psphot: $(BIN)/psphot.$(ARCH)

trunk/psphot/src/apply_psf_model.c

@@ -25 +25 @@
 
     float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
-    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    // float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
 
     // set the object surface-brightness limit for fitted pixels
@@ -66 +66 @@
         Nfit ++;
 
-        mark_psf_source (source);
+        mark_psf_source (source, shapeNsigma);
         subtract_psf_source (source);
     }

trunk/psphot/src/by_SN.c

@@ -1 +1 @@
 # include "psphot.h"
 
+// sort by SN (descending)
 int by_SN (const void **a, const void **b)
 {
-    psSource *A = *a;
-    psSource *B = *b;
+    psSource *A = *(psSource **)a;
+    psSource *B = *(psSource **)b;
 
     psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
@@ -10 +11 @@
 
     psF32 diff = fA - fB;
-    if (diff < FLT_EPSILON) return (-1);
-    if (diff > FLT_EPSILON) return (+1);
+    if (diff > FLT_EPSILON) return (-1);
+    if (diff < FLT_EPSILON) return (+1);
     return (0);
 }

trunk/psphot/src/choose_psf_model.c

@@ -20 +20 @@
         psArrayAdd (stars, 200, source);
     }
+    psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
 
     // define model fit pixels

trunk/psphot/src/mark_psf_source.c

@@ -14 +14 @@
 // any object which meets the criterion is marked as 
 // PS_SOURCE_BRIGHTSTAR 
-bool mark_psf_source (psArray *sources)
+bool mark_psf_source (psSource *source, float shapeNsigma)
 { 
-    bool  status      = false;
     float dSX, dSY, SX, SY, SN;
     float nSx, nSy;

trunk/psphot/src/psPolynomials.c

@@ -233 +233 @@
     // Build the B and A data structs.
     for (int k = 0; k < x->n; k++) {
-        if ((mask != NULL) && mask->U8[k]) continue;
+        if ((mask != NULL) && mask->data.U8[k]) continue;
         xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
       
@@ -325 +325 @@
     // Build the B and A data structs.
     for (int k  = 0; k < x->n; k++) {
-      if ((mask != NULL) && mask->data.U8[i]) continue;
+      if ((mask != NULL) && mask->data.U8[k]) continue;
         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
       

trunk/psphot/src/psUtils.c

@@ -261 +261 @@
 }
     
+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);
+}           

trunk/psphot/src/psphot.c

@@ -1 +1 @@
 # include "psphot.h"
+void dump_image (psImage *image, char *filename);
 
 int main (int argc, char **argv) {
@@ -33 +34 @@
     psf = choose_psf_model (image, config, sources);
 
+    psImage *zap = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
+    psImageInit (zap, 0);
+    find_defects (zap, sources, config, psf);
+    dump_image (zap, argv[2]);
+    exit (0);
+
     // test PSF on all except SATURATE and DEFECT (artifacts)
     apply_psf_model (image, config, sources, psf, sky);
@@ -44 +51 @@
     // subtract_galaxies (sources, config);
  
-    // write out subtracted image
-    unlink (argv[2]);
-    psFits *fits = psFitsAlloc (argv[2]);
-    psFitsWriteImage (fits, NULL, image, 0, NULL);
-    psFree (fits);
-
+    dump_image (image, argv[2]);
     DumpModelPSF (sources, "sources.dat");
     exit (0);
@@ -59 +61 @@
     exit (2);
 }
+
+void dump_image (psImage *image, char *filename) {
+
+    unlink (filename);
+    psFits *fits = psFitsAlloc (filename);
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFree (fits);
+    return;
+}
+

trunk/psphot/src/psphot.h

@@ -21 +21 @@
 psF64 psTimerMark (char *name);
 
+bool psImageInit (psImage *image,...);
 void psImageSmooth (psImage *image, float sigma, float Nsigma);
 psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
@@ -40 +41 @@
 psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
 psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
-psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, const psVector* x, const psVector* y, const psVector* yErr);
-psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
+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(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+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);
@@ -58 +60 @@
 pmPSF_Test *pmPSF_TestModel (psArray *stars, char *modelName);
 
-bool pmPSFFromModels (pmPSF *psf, psArray *models);
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
 psModel *psModelFromPSF (psModel *model, pmPSF *psf);
 void pmSourceMaskRegion (psSource *source, psRegion *region);
@@ -75 +77 @@
 bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
 bool subtract_galaxies (psArray *sources, psMetadata *config);
+
+int by_SN (const void **a, const void **b);
+bool subtract_psf_source (psSource *source);
+bool mark_psf_source (psSource *source, float shapeNsigma);
+bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);

trunk/psphot/src/pspsf.c

@@ -68 +68 @@
     }
     psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
     DumpModelFits (test->modelFLT, "modelsFLT.dat");
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
+    
+    // count valid sources
+    int Nkeep = 0;
+    for (int i = 0; i < sources->n; i++) {
+      if (test->mask->data.U8[i]) continue;
+      Nkeep++;
+    }
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
 
     // stage 3: refit with fixed shape parameters
@@ -96 +105 @@
     }
     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
     DumpModelFits (test->modelPSF, "modelsPSF.dat");
 
@@ -152 +162 @@
 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
 
-    int n;
-
     // construct the fit vectors from the collection of objects
     // use the mask to ignore missing fits 
@@ -187 +195 @@
         psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
         // psPolynomial2DDump (psf->params->data[i]);
+
+        // count valid sources
+        int Nkeep = 0;
+        for (int j = 0; j < mask->n; j++) {
+          if (mask->data.U8[j]) continue;
+          Nkeep++;
+        }
+        psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
+
     }
     return (true);

trunk/psphot/src/source_moments.c

@@ -18 +18 @@
     psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
     keep           = psRegionForImage (keep, image, keep);
-    peaks          = pmPeaksSubset (allpeaks, 100000, keep);
-
     // I need to exclude the saturated objects here, but include them again later...?
     peaks          = pmPeaksSubset (allpeaks, 100000, keep);

trunk/psphot/src/subtract_galaxies.c

@@ -5 +5 @@
 bool subtract_galaxies (psArray *sources, psMetadata *config) 
 { 
-    bool  status;
-
     for (int i = 0; i < sources->n; i++) {
         psSource *source = sources->data[i];

trunk/psphot/src/subtract_psf_source.c

@@ -4 +4 @@
 // need to control application of this with some thresholding
 
-bool subtract_psf_sources (psSource *source)
+bool subtract_psf_source (psSource *source)
 { 
   float sky;
-  psSource *source = sources->data[i];
 
   // non-stellar sources are ignored
@@ -15 +14 @@
   psModel *model = source->modelPSF;
   if (model == NULL) {
-    psLogMsg ("psphot.subtract_psf_source", "missing model for %f, %f\n", source->moments->x, source->moments->y);
+    psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y);
     return (false); 
   }         

trunk/psphot/src/test_psf_scatter.c

@@ -70 +70 @@
     snBin ->n = bin;
 
+    // XXX these APIs should get changed
     psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsxLine, NULL, snBin, dsxBin, NULL);
+    psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL);
     psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
-    psVectorFitPolynomial1D (dsyLine, NULL, snBin, dsyBin, NULL);
+    psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL);
     // these should have a slope of 1.0
     psPolynomial1DDump (dsxLine);