Changeset 37811

Timestamp:

Jan 11, 2015, 2:57:41 PM (11 years ago)

Author:

eugene

Message:

merge changes from ipp-20141224

Location:

trunk

Files:

• psLib/src/imageops/psImageMapFit.c (modified) (1 diff)
• psLib/src/math/psMinimizePolyFit.c (modified) (1 prop)
• psModules/src/objects/pmPSFtryMakePSF.c (modified) (2 diffs)
• psModules/src/objects/pmSourceIO_CMF.c.in (modified) (1 prop)
• psModules/src/objects/pmTrend2D.c (modified) (2 diffs)
• psphot/src/psphotChoosePSF.c (modified) (2 diffs)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (1 diff)
• psphot/src/psphotGalaxyParams.c (modified) (14 diffs)

trunk/psLib/src/imageops/psImageMapFit.c

@@ -286 +286 @@
     // fprintf (stderr, "Total: %f\n", Total);
 
+    double MaxPivot = 0.0;
+    for (int i = 0; i < Nx*Ny; i++) {
+      MaxPivot = PS_MAX(MaxPivot, fabs(A->data.F32[i][i]));
+      // fprintf (stderr, "piv, max: %f : %f\n", A->data.F32[i][i], MaxPivot);
+    }
+
     // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0
     psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8);
     psVectorInit (Empty, 0);
+    double MinPivot = 0.025*MaxPivot;
     for (int i = 0; i < Nx*Ny; i++) {
-        if (A->data.F32[i][i] == 0.0) {
+      if (fabs(A->data.F32[i][i]) < MinPivot) {
             Empty->data.S8[i] = 1;
             for (int j = 0; j < Nx*Ny; j++) {

trunk/psModules/src/objects/pmPSFtryMakePSF.c

@@ -45 +45 @@
 #include "pmSourceVisual.h"
 
+bool pmPSF_DataDump (char *filename, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *mask);
+
 /*****************************************************************************
 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of source->modelEXT entries
@@ -356 +358 @@
 }
 
+bool pmPSF_DataDump (char *filename, psVector *x, psVector *y, psVector *e0, psVector *e1, psVector *e2, psVector *mask) {
+
+
+  FILE *f = fopen (filename, "w");
+  if (!f) return false;
+
+  for (int i = 0; i < x->n; i++) {
+
+    fprintf (f, "%6.1f %6.1f : %5.2f %5.2f %5.2f : %2d\n", 
+             x->data.F32[i], 
+             y->data.F32[i], 
+             e0->data.F32[i], 
+             e1->data.F32[i], 
+             e2->data.F32[i], 
+             mask->data.U8[i]);
+  }
+  fclose (f);
+  return true;
+}

trunk/psModules/src/objects/pmTrend2D.c

@@ -211 +211 @@
         // XXX need to add the API which adjusts the scale
         status = psImageMapClipFit(pGoodFit, trend->map, trend->stats, mask, maskVal, x, y, f, df);
+        if (!status) {
+          psError(PS_ERR_UNKNOWN, true, "failed to build PSF model map");
+          return false;
+        }
+
+        // the psf model map can have nan pixels: repair this by extrapolation / interpolation
+        // XXX TEST: p_psImagePrint(0, trend->map->map, "before");
+        status = psImageMapRepair (trend->map->map);
+        // XXX TEST: p_psImagePrint(0, trend->map->map, "after");
+        if (!status) {
+          psError(PS_ERR_UNKNOWN, true, "failed to repair PSF model map");
+          return false;
+        }
         break;
 
@@ -217 +230 @@
     }
     return status;
+}
+
+bool psImageMapRepair (psImage *map) {
+
+
+  // XXX why is my repair not working??
+
+    // patch over bad regions (use average of 8 possible neighbor pixels)
+    // XXX consider testing all pixels against the 8 neighbors and replacing outliers...
+    double Count = 0;                   // number of good pixels
+    double Value = 0;                   // sum of good pixel's value
+    for (int iy = 0; iy < map->numRows; iy++) {
+        for (int ix = 0; ix < map->numCols; ix++) {
+            if (isfinite(map->data.F32[iy][ix])) {
+                Value += map->data.F32[iy][ix];
+                Count++;
+                continue;
+            }
+
+            double value = 0;
+            double count = 0;
+            for (int jy = iy - 1; jy <= iy + 1; jy++) {
+                if (jy <   0) continue;
+                if (jy >= map->numRows) continue;
+                for (int jx = ix - 1; jx <= ix + 1; jx++) {
+                    if (!jx && !jy) continue;
+                    if (jx   <   0) continue;
+                    if (jx   >= map->numCols) continue;
+                    if (!isfinite(map->data.F32[jy][jx])) continue;
+                    value += map->data.F32[jy][jx];
+                    count += 1.0;
+                }
+            }
+            if (count > 0) {
+              // psLogMsg ("psphot", PS_LOG_DETAIL, "patching image map %d, %d: %f (%d pts)\n", ix, iy, (value / count), (int) count);
+                map->data.F32[iy][ix] = value / count;
+            }
+        }
+    }
+    if (Count == 0) {
+        psError(PS_ERR_UNKNOWN, true, "failed to repair PSF model map");
+        return false;
+    }
+    Value /= Count;
+
+    // patch over remaining bad regions (use global average)
+    for (int iy = 0; iy < map->numRows; iy++) {
+        for (int ix = 0; ix < map->numCols; ix++) {
+            if (!isnan(map->data.F32[iy][ix])) continue;
+            map->data.F32[iy][ix] = Value;
+        }
+    }
+    return true;
 }
 

trunk/psphot/src/psphotChoosePSF.c

@@ -476 +476 @@
             psFree (modelPSF);
 
-            float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
-            fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
+            // float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
+            // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
         }
     }
@@ -630 +630 @@
         psFree (modelPSF);
 
-        float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
-        fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
+        // float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
+        // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
     }
 

trunk/psphot/src/psphotExtendedSourceAnalysis.c

@@ -38 +38 @@
 
 /*** for the moment, this test code : it is not thread safe ***/
-int    Nall = 0;
-int  Nskip1 = 0;
-int  Nskip2 = 0;
-int  Nskip3 = 0;
-int  Nskip4 = 0;
-int  Nskip5 = 0;
-int  Nskip6 = 0;
+static int    Nall = 0;
+static int  Nskip1 = 0;
+static int  Nskip2 = 0;
+static int  Nskip3 = 0;
+static int  Nskip4 = 0;
+static int  Nskip5 = 0;
+static int  Nskip6 = 0;
 
 # define SKIP(VALUE) { VALUE++; continue; }

trunk/psphot/src/psphotGalaxyParams.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-#define USE_FOOTPRINTS 1
+# define PRINT_STUFF 0
+# define USE_FOOTPRINTS 1
 
 // measure some properties on the exended objects
@@ -40 +41 @@
 }
 
-#ifdef notdef
+// # define COUNT_SKIPS
+# ifdef COUNT_SKIPS
 /*** for the moment, this test code : it is not thread safe ***/
-int  Nall = 0;
-int  Nskip1 = 0;
-int  Nskip2 = 0;
-int  Nskip3 = 0;
-int  Nskip4 = 0;
-int  Nskip5 = 0;
-int  Nskip6 = 0;
-
+static int  Nskip1 = 0;
+static int  Nskip2 = 0;
+static int  Nskip3 = 0;
+static int  Nskip4 = 0;
+static int  Nskip5 = 0;
+static int  Nskip6 = 0;
+static int  Nskip7 = 0;
+static int  Nskip8 = 0;
 # define SKIP(VALUE) { VALUE++; continue; }
-
-#endif //notdef
-
-#define SKIP(VALUE) {continue;}
+# else
+# define SKIP(VALUE) {continue;}
+# endif //notdef
 
 // aperture-like measurements for extended sources
@@ -104 +105 @@
 
     psImage *footprintImage = NULL;
-#ifdef USE_FOOTPRINTS
+#if USE_FOOTPRINTS
     footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     psImageInit(footprintImage, 0);
-    pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, true);
+    pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, false);
 #endif
 
@@ -207 +208 @@
 
 // # if (PS_TRACE_ON)
-# if (0)
+# ifdef COUNT_SKIPS
     fprintf (stderr, "ext analysis skipped @ 1  : %d\n", Nskip1);
     fprintf (stderr, "ext analysis skipped @ 2  : %d\n", Nskip2);
@@ -214 +215 @@
     fprintf (stderr, "ext analysis skipped @ 5  : %d\n", Nskip5);
     fprintf (stderr, "ext analysis skipped @ 6  : %d\n", Nskip6);
+    fprintf (stderr, "ext analysis skipped @ 7  : %d\n", Nskip7);
+    fprintf (stderr, "ext analysis skipped @ 8  : %d\n", Nskip8);
 #endif
 
@@ -221 +224 @@
 
     // psphotVisualShowPetrosians (sources);
+
+    // XXX TEST: SAVE IMAGE:
+    // psphotSaveImage (NULL, readout->image, "resid.fits");
 
     return true;
@@ -240 +246 @@
     psMetadata *recipe      = job->args->data[3];
 
-#ifndef USE_FOOTPRINTS
+#if (!USE_FOOTPRINTS)
     float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
 #endif
@@ -256 +262 @@
     // don't try and use bad models
     pmModelStatus badModel = PM_MODEL_STATUS_NONE;
-    badModel |= PM_MODEL_STATUS_NONCONVERGE;
+//    badModel |= PM_MODEL_STATUS_NONCONVERGE;  this causes most objects to get unmeasured
     badModel |= PM_MODEL_STATUS_BADARGS;
     badModel |= PM_MODEL_STATUS_OFFIMAGE;
@@ -291 +297 @@
 
         // Nall ++;
+
+        bool testObject = FALSE;
+        // testObject = testObject || ((fabs(source->peak->xf - 5788) < 10) && (fabs(source->peak->yf - 5711) < 10));
+        // testObject = testObject || ((fabs(source->peak->xf -   96) < 10) && (fabs(source->peak->yf - 5550) < 10));
+        // testObject = testObject || ((fabs(source->peak->xf -  300) < 10) && (fabs(source->peak->yf - 5988) < 10));
+        if (testObject) {
+          fprintf (stderr, "test object: %f, %f\n", source->peak->xf, source->peak->yf);
+          testObject = TRUE;
+        }
 
         // rules for measuring petrosian parameters for specific objects are set in
@@ -403 +418 @@
 #endif
 
-#ifdef USE_FOOTPRINTS
+#if (USE_FOOTPRINTS)
         psS32 **vFootprint = footprintImage->data.S32;
 #else
@@ -436 +451 @@
                 }
 
-#ifdef USE_FOOTPRINTS
-                if (vFootprint[yImage][xImage]) {
-#else
-                if (fabs(pixel) >= threshold) {
-#endif
+#if (USE_FOOTPRINTS)
+                // this is a pixel in the object
+                if (vFootprint[yImage][xImage] == source->peak->footprint->id) {
                     SET_MARK(MARK_SOURCE);
                     continue;
+                }
+                // this is a pixel in another object
+                if (vFootprint[yImage][xImage] > 0.0) continue;
+#else
+                // this is a pixel in an object (this or the other)
+                if (fabs(pixel) >= threshold) {
+                  continue
                 }
+#endif
                 // we have a background pixel.
                 // Check whether it's mirror is a background pixel as well.
@@ -456 +477 @@
                 if (isnan(mirrorPixel)) continue;
 
-#ifdef USE_FOOTPRINTS
+#if (USE_FOOTPRINTS)
                 if (vFootprint[yFlip + row0][xFlip + col0] == 0) {
 #else
@@ -581 +602 @@
         source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA;
 
-#ifdef PRINT_STUFF
+#if (PRINT_STUFF)
         if (Ngood % 40 == 1) {
-            printf("   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
-        }
-        printf("%5d %5.2e %3.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
-            source->id, sumI, source->extpars->ghalfLightRadius, radius, source->extpars->gRT, source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n);
+          fprintf (stderr, "   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
+        }
+        fprintf (stderr, "%5d %5.1f %5.1f : %5.2e %4.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
+                 source->id, source->peak->xf, source->peak->yf, 
+                 sumI, source->extpars->ghalfLightRadius, 
+                 radius, 
+                 source->extpars->gRT, 
+                 source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n);
 #endif
 
@@ -592 +617 @@
 restoreModel:
 
-        if ((extModel != sersicModel) || !isExtended) {
+        if (testObject) {
+          fprintf (stderr, "done with %d\n", source->id);
+        }
+
+# define USE_BEST_FIT TRUE
+        // XXX TEST: Keep the sersic model regardless
+        if (USE_BEST_FIT && ((extModel != sersicModel) || !isExtended)) {
             // add the sersic model back in
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {