Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c	(revision 34376)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c	(revision 34377)
@@ -10,10 +10,6 @@
 
 bool psImageQuickStats (QuickStats *stats, psImage *image, psRegion region);
-
-bool  psphotTestRadialModel (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logFlux, psVector *logRad, float Rmax);
-bool  psphotTestRadialModelSub (pmSource *source, psVector *logRmodel, psVector *logFmodel, float Xo, float Yo, float Rmax, psImageMaskType maskVal);
-bool  psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal);
 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
-bool  psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel);
+bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -46,7 +42,9 @@
     TBD:
 
-    1) carry the radial profile out of this function (element of pmSource)
-    2) function to replace the radial profile for a source
-    3) recenter the profile (based)
+    * function to replace the radial profile for a source
+    * recenter the profile (based on cross-correlation / convolution with 1D profile)
+
+    * raise a bit somewhere for these super saturated stars (if they were so modeled)
+    * consider the subtraction bit : when to raise it?
 
  **/
@@ -171,8 +169,13 @@
 	source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
 
+	if (!psphotSatstarProfileModel (source, maskVal)) continue;
+
+	source->mode |= PM_SOURCE_MODE_SATSTAR; // yes, this source IS saturated
+	source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile
+
 	// XXX visualize, model, and subtract
-	if (!psphotVisualRadialProfileSatstar (source, maskVal)) {
-	  break;
-	}
+	// if (!psphotVisualRadialProfileSatstar (source, maskVal)) {
+	// break;
+	// }
 
 	// generate radial profile, store on the source structure
@@ -392,4 +395,252 @@
 }
 
+// create a model for the radial profile of a saturated star (is this actually more generic?)
+bool psphotSatstarProfileModel (pmSource *source, psImageMaskType maskVal) {
+
+    // XXX user define somewhere?
+    float Rmax = 320.0;
+
+    // XXX is this ever the case??
+    bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
+    if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+
+    int nPts = source->pixels->numRows * source->pixels->numCols;
+    psVector *logR = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *flux = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+
+    int ng = 0;
+
+    // choose the best center for this profile
+    float Xo = NAN;
+    float Yo = NAN;
+
+    if (source->modelPSF) {
+	// XXX do we ever have a PSF model for a SATSTAR?
+	Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+	Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+	Xo = source->moments->Mx;
+	Yo = source->moments->My;
+    }
+
+    float Xc = Xo - source->pixels->col0 - 0.5;
+    float Yc = Yo - source->pixels->row0 - 0.5;
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+	for (int ix = 0; ix < source->pixels->numCols; ix++) {
+	    if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue;
+	    // XXX do this faster by generating R^2 and returning 0.5*log10(R^2)?
+	    logR->data.F32[ng] = log10(hypot (ix - Xc, iy - Yc));
+	    flux->data.F32[ng] = source->pixels->data.F32[iy][ix];
+	    ng++;
+	}
+    }
+    logR->n = ng;
+    flux->n = ng;
+
+    // XXX do something sensible here if there are no pixels
+
+    // measure the radial profile 
+    psVector *logFmodel = NULL;
+    psVector *logRmodel = NULL;
+    psphotSatstarProfileCreate (source, &logRmodel, &logFmodel, logR, flux, Rmax);
+    
+    // XXX do something sensible here if the profile is crap
+
+    source->satstar = pmSourceSatstarAlloc();
+    source->satstar->Xo = Xo;
+    source->satstar->Yo = Yo;
+    source->satstar->Rmax = Rmax;
+    source->satstar->logFmodel = logFmodel;
+    source->satstar->logRmodel = logRmodel;
+
+    // subtract the profile (false => subtract)
+    psphotSatstarProfileOp (source, maskVal, 1.0, 0, false);
+
+    psFree (logR);
+    psFree (flux);
+
+    return true;
+}
+
+// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut
+bool psphotSatstarProfileCreate (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logR, psVector *flux, float Rmax) {
+
+  // we have log(radius) & log(flux).  find the median flux in log radial bins
+
+  float logRmax = log10(Rmax);
+  float logRdel = 0.1; // XXX user-define?
+
+  psVector *logRmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
+  psVector *logFmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
+
+  pmSourceRadialProfileSortPair (logR, flux);
+
+  psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+  psVector *fluxVals = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+  int bin = 0;
+  for (int i = 0; i < logRmodel->n; i++) {
+    
+    // bin (i) has log radius range (i*logRdel : (i+1)*logRdel)
+    float lRmin = logRdel*(i + 0);
+    float lRmax = logRdel*(i + 1);
+
+    // reset the flux vector
+    fluxVals->n = 0;
+    psStatsInit (fluxStats);
+
+    while (logR->data.F32[bin] < lRmax) {
+      if (isfinite(flux->data.F32[bin])) {
+	  psVectorAppend (fluxVals, flux->data.F32[bin]);
+      }
+      bin ++;
+    }
+    
+    // we have the set of fluxes for this bin; find the median values
+    
+    float Rmin = pow(10.0, lRmin);
+    float Rmax = pow(10.0, lRmax);
+
+    float Rmean = (2.0/3.0) * (pow(Rmax, 3.0) - pow(Rmin, 3.0)) / (PS_SQR(Rmax) - PS_SQR(Rmin));
+    logRmodel->data.F32[i] = log10(Rmean);
+
+    float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin));
+    if (fluxVals->n < 0.25*Area) {
+      logFmodel->data.F32[i] = NAN;
+      continue;
+    }
+    
+    psVectorStats (fluxStats, fluxVals, NULL, NULL, 0);
+    if (fluxStats->robustMedian > 0.0) {
+	logFmodel->data.F32[i] = log10(fluxStats->robustMedian);
+    } else {
+	logFmodel->data.F32[i] = -3.0;
+    }
+    // fprintf (stderr, "R: %f, F: %f +/- %f\n", Rmean, fluxStats->robustMedian, fluxStats->robustStdev);
+  }
+
+  // now how do i use this to subtract a model??
+  *logRmodelOut = logRmodel;
+  *logFmodelOut = logFmodel;
+  
+  // need to free stuff
+  psFree (fluxStats);
+  psFree (fluxVals);
+
+  return true;
+}
+
+bool psphotSatstarProfileOp (pmSource *source, psImageMaskType maskVal, float FACTOR, pmModelOpMode mode, bool add) {
+
+    int alt;
+    float logRdel = 0.1;
+
+    float Xc = source->satstar->Xo - source->pixels->col0 - 0.5;
+    float Yc = source->satstar->Yo - source->pixels->row0 - 0.5;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+	for (int ix = 0; ix < source->pixels->numCols; ix++) {
+	    if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue;
+
+	    float radius = hypot (ix - Xc, iy - Yc) ;
+	    float logR = log10(radius);
+
+	    int bin = (int)(logR / logRdel);
+
+	    // we are going to interpolate if possible, or extrapolate if necessary
+	    if (bin >= logRmodel->n - 1) bin = logRmodel->n - 1;
+	    if (bin < 0) bin = 0;
+      
+	    // interpolate to the current radial position
+	    // XXX BIG HACK : skip nan bins for now
+
+	    float logF0 = logFmodel->data.F32[bin];
+	    float logR0 = logRmodel->data.F32[bin];
+	    if (!isfinite(logF0)) continue; 
+    
+	    // interpolate between closest two bins if possible, extrapolate on ends
+	    if (logR < logR0) {
+		alt = (bin > 0) ? bin - 1 : bin + 1;
+	    } else {
+		alt = (bin < logRmodel->n - 1) ? bin + 1 : bin - 1;
+	    }
+
+	    float logF1 = logFmodel->data.F32[alt];
+	    float logR1 = logRmodel->data.F32[alt];
+	    if (!isfinite(logF1)) continue; 
+
+	    // XXX use linear flux, not logFlux
+	    float logF = InterpolateValues (logR0, logF0, logR1, logF1, logR);
+	    float flux = pow (10.0, logF);
+
+	    if (mode & PM_MODEL_OP_NOISE) {
+		if (add) {
+		    source->variance->data.F32[iy][ix] += FACTOR * flux;
+		} else {
+		    source->variance->data.F32[iy][ix] -= FACTOR * flux;
+		}
+	    } else {
+		if (add) {
+		    source->pixels->data.F32[iy][ix] += flux;
+		} else {
+		    source->pixels->data.F32[iy][ix] -= flux;
+		}
+	    }
+	}
+    }
+    return true;
+}
+
+// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut
+bool psphotSatstarPhotometry (pmSource *source) {
+
+    if (!source->satstar) return false;
+
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+    
+    float fluxTotal = 0.0;
+    float logRdel = 0.1; // XXX user-define (or carry in satstar)
+
+    // integrate flux in radial bins 
+    for (int i = 0; i < logRmodel->n; i++) {
+	// just add up the mean flux in each bin and multiply by the area of the bin
+	
+	float logF = logFmodel->data.F32[i];
+	if (!isfinite(logF)) continue;
+	float flux = pow(10.0, logF); // this is the mean flux per pixel (ie, surface brightness)
+
+	// bin (i) has log radius range (i*logRdel : (i+1)*logRdel)
+	float lRmin = logRdel*(i + 0);
+	float lRmax = logRdel*(i + 1);
+
+	float Rmin = pow(10.0, lRmin);
+	float Rmax = pow(10.0, lRmax);
+
+	float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin));
+
+	fluxTotal += flux * Area;
+    }
+
+    source->psfMag    = -2.5*log10(fluxTotal);
+    source->psfMagErr = 0.05; 
+
+    // XXX I have no idea of a realistic error on the photometry here.  the bottom line is that
+    // the error is totally dominated by the loss of charge due to saturation.  I am not
+    // modeling the profile in any real detail other than to follow the radial bins.  
+
+    source->extMag    = NAN;
+    source->apMag     = NAN;
+    source->apMagRaw  = NAN;
+    source->apFlux    = fluxTotal;
+    source->apFluxErr = 0.05*fluxTotal;
+
+    return true;
+}
+
+# if (0)
 static bool skipDisplay = false;
 bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal) {
@@ -743,118 +994,5 @@
     return true;
 }
-
-bool psphotTestRadialModel (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logRad, psVector *logFlux, float Rmax) {
-
-  // we have log(radius) & log(flux).  find the median flux in log radial bins
-
-  float logRmax = log10(Rmax);
-  float logRdel = 0.1;
-
-  psVector *logRmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
-  psVector *logFmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
-
-  pmSourceRadialProfileSortPair (logRad, logFlux);
-
-  psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-  psVector *flux = psVectorAllocEmpty (100, PS_TYPE_F32);
-
-  int bin = 0;
-  for (int i = 0; i < logRmodel->n; i++) {
-    
-    // bin (i) has log radius range (i*logRdel : (i+1)*logRdel)
-    float lRmin = logRdel*(i + 0);
-    float lRmax = logRdel*(i + 1);
-
-    // reset the flux vector
-    flux->n = 0;
-    psStatsInit (fluxStats);
-
-    while (logRad->data.F32[bin] < lRmax) {
-      if (isfinite(logFlux->data.F32[bin])) {
-	  psVectorAppend (flux, logFlux->data.F32[bin]);
-      }
-      bin ++;
-    }
-    
-    // we have the set of fluxes for this bin; find the median values
-    
-    float Rmin = pow(10.0, lRmin);
-    float Rmax = pow(10.0, lRmax);
-
-    float Rmean = (2.0/3.0) * (pow(Rmax, 3.0) - pow(Rmin, 3.0)) / (PS_SQR(Rmax) - PS_SQR(Rmin));
-    logRmodel->data.F32[i] = log10(Rmean);
-
-    float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin));
-    if (flux->n < 0.25*Area) {
-      logFmodel->data.F32[i] = NAN;
-      continue;
-    }
-    
-    psVectorStats (fluxStats, flux, NULL, NULL, 0);
-    if (fluxStats->robustMedian > 0.0) {
-	logFmodel->data.F32[i] = log10(fluxStats->robustMedian);
-    } else {
-	logFmodel->data.F32[i] = -3.0;
-    }
-    // fprintf (stderr, "R: %f, F: %f +/- %f\n", Rmean, fluxStats->robustMedian, fluxStats->robustStdev);
-  }
-
-  // now how do i use this to subtract a model??
-  *logRmodelOut = logRmodel;
-  *logFmodelOut = logFmodel;
-  
-  // need to free stuff
-  psFree (fluxStats);
-  psFree (flux);
-
-  return true;
-}
-
-bool psphotTestRadialModelSub (pmSource *source, psVector *logRmodel, psVector *logFmodel, float Xo, float Yo, float Rmax, psImageMaskType maskVal) {
-
-    int alt;
-  float logRdel = 0.1;
-
-  for (int iy = 0; iy < source->pixels->numRows; iy++) {
-    for (int ix = 0; ix < source->pixels->numCols; ix++) {
-      if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue;
-
-      float radius = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
-      float logR = log10(radius);
-
-      int bin = (int)(logR / logRdel);
-
-      // we are going to interpolate if possible, or extrapolate if necessary
-      if (bin >= logRmodel->n - 1) bin = logRmodel->n - 1;
-      if (bin < 0) bin = 0;
-      
-      // interpolate to the current radial position
-      // XXX BIG HACK : skip nan bins for now
-
-      float logF0 = logFmodel->data.F32[bin];
-      float logR0 = logRmodel->data.F32[bin];
-      if (!isfinite(logF0)) continue; 
-    
-      // interpolate between closest two bins if possible, extrapolate on ends
-      if (logR < logR0) {
-	  alt = (bin > 0) ? bin - 1 : bin + 1;
-      } else {
-	  alt = (bin < logRmodel->n - 1) ? bin + 1 : bin - 1;
-      }
-
-      float logF1 = logFmodel->data.F32[alt];
-      float logR1 = logRmodel->data.F32[alt];
-      if (!isfinite(logF1)) continue; 
-
-      // XXX use linear flux, not logFlux
-      float logF = InterpolateValues (logR0, logF0, logR1, logF1, logR);
-      float flux = pow (10.0, logF);
-
-      // float flux = InterpolateValues (logR0, logF0, logR1, logF1, logR);
-      source->pixels->data.F32[iy][ix] -= flux;
-    }
-  }
-  return true;
-}
+# endif
 
 // only valid for F32
@@ -902,2 +1040,44 @@
     return true;
 }
+
+bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add) {
+
+    bool status;
+
+    psTimerStart ("psphot.deblend.sat");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
+	return true;
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // examine sources in decreasing SN order
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+	if (!(source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE)) continue;
+
+	if (!psphotSatstarProfileOp (source, maskVal, 1.0, 0, add)) continue;
+    }
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "satstar op: %f sec\n", psTimerMark ("psphot.deblend.sat"));
+    return true;
+}
+