Changeset 25128

Timestamp:

Aug 19, 2009, 12:16:54 PM (17 years ago)

Author:

eugene

Message:

updates to make radial profiles a bit more robust

Location:

branches/eam_branches/20090715/psphot/src

Files:

• psphotPetrosianProfile.c (modified) (1 diff)
• psphotPetrosianStudy.c (modified) (6 diffs)
• psphotPetrosianVisual.c (modified) (2 diffs)
• psphotRadialProfileByAngle.c (modified) (3 diffs)
• psphotRadiiFromProfiles.c (modified) (3 diffs)

branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c

@@ -13 +13 @@
 
   int Nsec = 24;
-  float Rmax = 64;
+  float Rmax = 200;
   float fluxMin = 0.0;
-  float fluxMax = 1000.0;
+  float fluxMax = source->peak->flux;
 
   // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux

branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+# define DX 512
+# define DY 512
+
+// XXX add noise and seeing.
+// XXX double check on sersic functional form
+// XXX modify ratio if ratio > 1.0 (swap major and minor)
+
+pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo);
 
 int main (int argc, char **argv) {
@@ -7 +16 @@
 
   int N;
-
-  // XXX add noise and seeing.
-  // XXX double check on sersic functional form
-  // XXX modify ratio if ratio > 1.0 (swap major and minor)
+  float Xo = 0.5*DX;
+  float Yo = 0.5*DY;
+  char *image = NULL;
+  pmSource *source = NULL;
 
   float peak = 1000.0;
@@ -18 +27 @@
   float sersic = 0.5;
   float skynoise = 0.0;
-
+  
   if ((N = psArgumentGet (argc, argv, "-peak"))) {
     psArgumentRemove (N, &argc, argv);
@@ -53 +62 @@
     pmVisualSetVisual(true);
   }
+  if ((N = psArgumentGet (argc, argv, "-coords"))) {
+    psArgumentRemove (N, &argc, argv);
+    Xo = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+    Yo = atof(argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
+  if ((N = psArgumentGet (argc, argv, "-image"))) {
+    psArgumentRemove (N, &argc, argv);
+    image = psStringCopy (argv[N]);
+    psArgumentRemove (N, &argc, argv);
+  }
 
   if (argc != 1) {
@@ -61 +82 @@
   // create a containing image & associated readout
   pmReadout *readout = pmReadoutAlloc(NULL);
-  readout->image = psImageAlloc(128, 128, PS_TYPE_F32);
+  if (!image) {
+      readout->image = psImageAlloc(DX, DY, PS_TYPE_F32);
 
-  // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0
-  readout->variance = psImageAlloc(128, 128, PS_TYPE_F32);
+      // create a dummy variance, but don't populate -- it is not used by pmSourceMoments if the sigma parameters is set to 0.0
+      readout->variance = psImageAlloc(DX, DY, PS_TYPE_F32);
 
-  // create a model & associated source
-  pmModelType type = pmModelClassGetType("PS_MODEL_SERSIC");
-  pmModel *model = pmModelAlloc(type);
+      // create a model & associated source
+      pmModelType type = pmModelClassGetType("PS_MODEL_SERSIC");
+      pmModel *model = pmModelAlloc(type);
 
-  // set the model parameters
-  model->params->data.F32[PM_PAR_SKY]  = 0.0;
-  model->params->data.F32[PM_PAR_I0]   = peak;
-  model->params->data.F32[PM_PAR_XPOS] = 64.0;
-  model->params->data.F32[PM_PAR_YPOS] = 64.0;
+      // set the model parameters
+      model->params->data.F32[PM_PAR_SKY]  = 0.0;
+      model->params->data.F32[PM_PAR_I0]   = peak;
+      model->params->data.F32[PM_PAR_XPOS] = Xo;
+      model->params->data.F32[PM_PAR_YPOS] = Yo;
 
-  psEllipseAxes axes;
-  axes.major = sigma;
-  axes.minor = sigma*ARatio;
-  axes.theta = angle;
+      psEllipseAxes axes;
+      axes.major = sigma;
+      axes.minor = sigma*ARatio;
+      axes.theta = angle;
 
-  psEllipseShape shape = psEllipseAxesToShape (axes);
+      psEllipseShape shape = psEllipseAxesToShape (axes);
 
-  // XXX set the sigma with user input
-  model->params->data.F32[PM_PAR_SXX]  = shape.sx * M_SQRT2;
-  model->params->data.F32[PM_PAR_SYY]  = shape.sy * M_SQRT2;
-  model->params->data.F32[PM_PAR_SXY]  = shape.sxy;
+      // XXX set the sigma with user input
+      model->params->data.F32[PM_PAR_SXX]  = shape.sx * M_SQRT2;
+      model->params->data.F32[PM_PAR_SYY]  = shape.sy * M_SQRT2;
+      model->params->data.F32[PM_PAR_SXY]  = shape.sxy;
 
-  if (model->params->n > 7) {
-    model->params->data.F32[PM_PAR_7]  = sersic;
+      if (model->params->n > 7) {
+          model->params->data.F32[PM_PAR_7]  = sersic;
+      }
+
+      // generate source container & populate image
+      source = pmSourceFromModel(model, readout, Xo, PM_SOURCE_TYPE_STAR);
+
+      // generate the modelFlux 
+      pmSourceCacheModel(source, 0);
+
+      // instantiate the source
+      pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); 
+
+      // XXX add noise here...
+      psphotSaveImage(NULL, readout->image, "sersic.fits");
+
+  } else {
+      psRegion full = psRegionSet(0,0,0,0);
+      psFits *fits = psFitsOpen(image, "r");
+      readout->image = psFitsReadImage(fits, full, 0);
+
+      source = pmSourceAlloc();
+      source->peak = psphotLocalPeak(readout, Xo, Yo);
+      pmSourceDefinePixels (source, readout, Xo, Yo, 64);
   }
-
-  pmSource *source = pmSourceFromModel(model, readout, 32.0, PM_SOURCE_TYPE_STAR);
-
-  // generate the modelFlux 
-  pmSourceCacheModel(source, 0);
-
-  // instantiate the source
-  pmSourceAdd(source, PM_MODEL_OP_FUNC, 0); 
-
-  // XXX add noise here...
-
-  psphotSaveImage(NULL, readout->image, "sersic.fits");
 
   psphotPetrosianProfile (source);
@@ -110 +142 @@
   exit (0);
 }
+
+// Xo, Yo are in parent coords
+pmPeak *psphotLocalPeak(pmReadout *readout, int Xo, int Yo) {
+
+    int Xp = Xo;
+    int Yp = Yo;
+    float peakFlux = readout->image->data.F32[Yp][Xp];
+
+    // find local peak within +/- 3 pix of the given coordinate
+    for (int iy = Yo - 3; iy <= Yo + 3; iy++) {
+        for (int ix = Xo - 3; ix <= Xo + 3; ix++) {
+            if (peakFlux < readout->image->data.F32[iy][ix]) {
+                Xp = ix;
+                Yp = iy;
+                peakFlux = readout->image->data.F32[Yp][Xp];
+            }
+        }
+    }
+
+    pmPeak *peak = pmPeakAlloc(Xp, Yp, peakFlux, PM_PEAK_LONE);
+
+    // calculate fractional peak position relative to Xp,Yp
+    psPolynomial2D *bicube = psImageBicubeFit (readout->image, Xp, Yp);
+    psPlane min = psImageBicubeMin (bicube);
+    psFree (bicube);
+
+    // if min point is too deviant, use the peak value
+    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
+        peak->xf = min.x + Xp;
+        peak->yf = min.y + Yp;
+
+        // xf,yf must land on image with 0 pixel border
+        peak->xf = PS_MAX (PS_MIN (peak->xf, readout->image->numCols - 1), readout->image->col0);
+        peak->yf = PS_MAX (PS_MIN (peak->yf, readout->image->numRows - 1), readout->image->row0);
+    } else {
+        peak->xf = Xp;
+        peak->yf = Yp;
+    }
+
+    return peak;
+}

branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c

@@ -26 +26 @@
         graphdata->xmin = graphdata->xmax = xVec->data.F32[0];
         for (int i = 1; i < xVec->n; i++) {
+            if (!isfinite(xVec->data.F32[i])) continue;
             graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
             graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
@@ -36 +37 @@
         graphdata->ymin = graphdata->ymax = yVec->data.F32[0];
         for (int i = 1; i < yVec->n; i++) {
+            if (!isfinite(yVec->data.F32[i])) continue;
             graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
             graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);

branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c

@@ -5 +5 @@
 
 // XXX choices for rMax and Nsec?
-// XXX where does the output go?
-// XXX how does this fit in context?
 
 bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
 
+    // we want to have an even number of sectors so we can do 180 deg symmetrizing
+    Nsec = (Nsec % 2) ? Nsec + 1 : Nsec;
     float dtheta = 2.0*M_PI / Nsec;
 
@@ -46 +46 @@
             // value is NAN if we run off the image
             float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
-            if (isnan(value)) continue;
+
+            // keep the nan values so all vectors are matched
+            // if (isnan(value)) continue;
 
             psVectorAppend (radius, r);
@@ -62 +64 @@
     }
 
+    for (int i = 0; i < Nsec / 2; i++) {
+
+        psVector *r1 = petrosian->radii->data[i];
+        psVector *r2 = petrosian->radii->data[i+Nsec/2];
+
+        psVector *f1 = petrosian->fluxes->data[i];
+        psVector *f2 = petrosian->fluxes->data[i+Nsec/2];
+
+        psAssert (r1->n == r2->n, "mis-matched vectors");
+        psAssert (f1->n == f2->n, "mis-matched vectors");
+
+        // we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair
+        for (int j = 0; j < r1->n; j++) {
+            
+            float flux;
+
+            if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) {
+                flux = NAN;
+                goto setflux;
+            }
+
+            if (!isfinite(f1->data.F32[j])) {
+                flux = f2->data.F32[j];
+                goto setflux;
+            }
+            if (!isfinite(f2->data.F32[j])) {
+                flux = f1->data.F32[j];
+                goto setflux;
+            }
+
+            flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]);
+
+        setflux:
+            f1->data.F32[j] = flux;
+            f2->data.F32[j] = flux;
+        }
+    }    
     return true;
 }

branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

@@ -30 +30 @@
 
     // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size
+    // XXX reconsider the fractional isophote value
     float Fm = fluxMin + 0.25*fluxRange;
     float Fp = fluxMin + 0.75*fluxRange;
     float Fo = fluxMin + 0.50*fluxRange;
+    int Rbin = 1;
       
-    // find the mean radius of the points in the flux range Fm - Fp:
-    float Rsum = 0;
-    float Rnpt = 0;
-    for (int i = 0; i < flux->n; i++) {
-        if (flux->data.F32[i] < Fm) continue;
-        if (flux->data.F32[i] > Fp) continue;
-      
-        Rsum += radius->data.F32[i];
-        Rnpt ++;
-    }
+    // find the median radius of the points in the flux range Fm - Fp:
+    { 
+        // storage vector for stats
+        psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
 
-    // rebin with bins 1/2 the size of the mean radius
-    int Rbin = 0;
-    if (Rnpt == 0) {
-        Rbin = 1;
-    } else {
-        Rbin = MAX(1, 0.5*(Rsum / Rnpt));
+        for (int i = 0; i < flux->n; i++) {
+            if (!isfinite(flux->data.F32[i])) continue;
+            if (flux->data.F32[i] < Fm) continue;
+            if (flux->data.F32[i] > Fp) continue;
+            
+            psVectorAppend (values, radius->data.F32[i]);
+        }
+        psVectorStats (stats, values, NULL, NULL, 0);
+
+        // if we have a valid range, rebin with bin size 1/2 of median radius
+        if (isfinite(stats->sampleMedian)) {
+            Rbin = MAX(1, 0.5*stats->sampleMedian);
+        }
     }
 
@@ -94 +98 @@
                 psStatsInit(stats);
             }
+            if (!isfinite(flux->data.F32[i])) continue;
             psVectorAppend (values, flux->data.F32[i]);
         }
@@ -105 +110 @@
     bool above = true;
     for (int i = 0; i < fluxBinned->n; i++) {
+
+        if (!isfinite(fluxBinned->data.F32[i])) continue;
 
         // find the largest radius that matches the flux transition