Changeset 15562 for trunk/psphot/src

Timestamp:

Nov 9, 2007, 3:09:20 PM (19 years ago)

Author:

eugene

Message:

updating all changes from from eam_branch_20071023

Location:

trunk/psphot/src

Files:

• Makefile.am (modified) (1 diff)
• psphot.h (modified) (1 diff)
• psphotAnnuli.c (modified) (1 diff)
• psphotExtendedSources.c (modified) (1 diff)
• psphotIsophotal.c (modified) (1 diff)
• psphotMakeResiduals.c (modified) (3 diffs)
• psphotPetrosian.c (modified) (1 diff)
• psphotRadialProfile.c (modified) (3 diffs)

trunk/psphot/src/Makefile.am

@@ -47 +47 @@
         psphotModelWithPSF.c     \
         psphotExtendedSources.c  \
+        psphotRadialProfile.c    \
         psphotPetrosian.c        \
         psphotIsophotal.c        \

trunk/psphot/src/psphot.h

@@ -116 +116 @@
 psKernel       *psphotKernelFromPSF (pmSource *source, int nPix);
 
+bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool            psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal);
 bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal);

trunk/psphot/src/psphotAnnuli.c

@@ -3 +3 @@
 bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  bool status;
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *weight = source->extpars->profile->weight;
+  psVector *flux = source->extpars->profile->flux;
+
+  // XXX how do I define the radii?  we can put a vector in the recipe...
+  // radialBins defines the bounds or start and stop (we can skip some that way...
+  psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+  psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+  assert (radialBinsLower->n == radialBinsUpper->n);
+
+  psVector *fluxValues = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *fluxSquare = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *fluxWeight = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVector *pixelCount = psVectorAlloc (radialBinsLower->n, PS_TYPE_F32);
+  psVectorInit (fluxValues, 0.0);
+  psVectorInit (fluxSquare, 0.0);
+  psVectorInit (fluxWeight, 0.0);
+  psVectorInit (pixelCount, 0.0);
+
+  // XXX this code assumes the radii are in pixels.  convert from arcsec with plate scale
+  // XXX assume the annulii above are not overlapping?  much faster...
+  // XXX this might be must faster in the reverse order: loop over annulii and use disection to 
+  // skip to the start of the annulus.
+  for (int i = 0; i < flux->n; i++) {
+    for (int j = 0; j < radialBinsLower->n; j++) {
+      if (radius->data.F32[i] < radialBinsLower->data.F32[j]) continue;
+      if (radius->data.F32[i] > radialBinsUpper->data.F32[j]) continue;
+      fluxValues->data.F32[j] += flux->data.F32[i];
+      fluxSquare->data.F32[j] += PS_SQR(flux->data.F32[i]);
+      fluxWeight->data.F32[j] += weight->data.F32[i];
+      pixelCount->data.F32[j] += 1.0;
+    }
+  }
+
+  for (int j = 0; j < radialBinsLower->n; j++) {
+    fluxValues->data.F32[j] /= pixelCount->data.F32[j];
+    fluxSquare->data.F32[j] /= pixelCount->data.F32[j];
+    fluxSquare->data.F32[j] -= PS_SQR(fluxValues->data.F32[j]);
+  }
+  
+  source->extpars->annuli = pmSourceAnnuliAlloc ();
+  source->extpars->annuli->flux    = fluxValues;
+  source->extpars->annuli->fluxErr = fluxWeight;
+  source->extpars->annuli->fluxVar = fluxSquare;
+
   return true;
+}
 
-}

trunk/psphot/src/psphotExtendedSources.c

@@ -66 +66 @@
             return false;
         } else {
+          // XXX why am I caching the model?
             pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
             psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
             Npsf ++;
+        }
+
+        // all of these below require the radial profile
+        // XXX push this as a test and call in each of the functions below?
+        // XXX this constructs a pmSourceExtendedParameters element
+        if (doPetrosian || doIsophotal || doAnnuli || doKron) {
+          if (!psphotRadialProfile (source, recipe, maskVal)) {
+            psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile");
+            return false;
+          }
         }
 

trunk/psphot/src/psphotIsophotal.c

@@ -3 +3 @@
 bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  bool status;
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *flux = source->extpars->profile->flux;
+
+  // flux at which to measure isophotal parameters
+  // XXX cache this?
+  float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX");
+
+  // find the first bin below the flux level and the last above the level
+  // XXX can this be done faster with disection?
+  // XXX do I need to worry about crazy outliers?  
+  // XXX should i be smoothing or fitting the curve?
+  int firstBelow = -1;
+  int lastAbove = -1;
+  for (int i = 0; i < flux->n; i++) {
+    if (flux->data.F32[i] > ISOPHOT_FLUX) lastAbove = i;
+    if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i;
+  }
+
+  // need to examine pixels in this vicinity
+  float isophotalFluxFirst = 0;
+  float isophotalFluxLast = 0;
+  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
+    if (i <= firstBelow) {
+      isophotalFluxFirst += flux->data.F32[i];
+    }
+    if (i <= lastAbove) {
+      isophotalFluxLast += flux->data.F32[i];
+    }
+  }
+  float isophotalFlux    = 0.5*(isophotalFluxLast + isophotalFluxFirst);
+  float isophotalFluxErr = 0.5*fabs(isophotalFluxLast - isophotalFluxFirst);
+
+  float isophotalRad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
+  float isophotalRadErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
+
+  if (!source->extpars->isophot) {
+    source->extpars->isophot = pmSourceIsophotalValuesAlloc ();
+  }
+  
+  // these are uncalibrated: instrumental mags and pixel units
+  source->extpars->isophot->mag    = -2.5*log10(isophotalFlux);
+  source->extpars->isophot->magErr = isophotalFluxErr / isophotalFlux;
+
+  source->extpars->isophot->rad    = isophotalRad;
+  source->extpars->isophot->radErr = isophotalRadErr;
+
   return true;
 

trunk/psphot/src/psphotMakeResiduals.c

@@ -28 +28 @@
 
     float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
+    PS_ASSERT (status, false);
+
+    float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
     PS_ASSERT (status, false);
 
@@ -194 +197 @@
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
                 //resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
+
+                if (resid->Ro->data.F32[oy][ox] < pixelSN*fluxStats->sampleStdev) {
+                  resid->mask->data.U8[oy][ox] = 1;
+                }
+
             } else {
                 assert (SPATIAL_ORDER == 1);
@@ -227 +235 @@
                 resid->Rx->data.F32[oy][ox] = B->data.F64[1];
                 resid->Ry->data.F32[oy][ox] = B->data.F64[2];
+
+                float dRo = sqrt(A->data.F32[0][0]);
+                if (resid->Ro->data.F32[oy][ox] < pixelSN*dRo) {
+                  resid->mask->data.U8[oy][ox] = 1;
+                }
                 //resid->weight->data.F32[oy][ox] = XXX;
             }

trunk/psphot/src/psphotPetrosian.c

@@ -3 +3 @@
 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
 
-  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  bool status;
+
+  assert (source->extpars);
+  assert (source->extpars->profile);
+  assert (source->extpars->profile->radius);
+  assert (source->extpars->profile->flux);
+
+  psVector *radius = source->extpars->profile->radius;
+  psVector *flux = source->extpars->profile->flux;
+
+  // flux at which to measure isophotal parameters
+  // XXX cache this?
+  float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0");
+  float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO");
+
+  // first find flux at R0
+  int firstBelow = -1;
+  int lastAbove = -1;
+  for (int i = 0; i < radius->n; i++) {
+    if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i;
+    if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i;
+  }
+
+  // average flux in this range
+  float fluxR0 = 0.0;
+  int fluxRn = 0;
+  for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {
+    fluxR0 += flux->data.F32[i];
+    fluxRn ++;
+  }
+  fluxR0 /= (float)(fluxRn);
+  
+  // target flux for petrosian radius
+  float fluxRP = fluxR0 * PETROSIAN_RF;
+
+  // find the first bin below the flux level and the last above the level
+  // XXX can this be done faster with disection?
+  // XXX do I need to worry about crazy outliers?  
+  // XXX should i be smoothing or fitting the curve?
+  firstBelow = -1;
+  lastAbove = -1;
+  for (int i = 0; i < flux->n; i++) {
+    if (flux->data.F32[i] > fluxRP) lastAbove = i;
+    if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i;
+  }
+
+  // need to examine pixels in this vicinity
+  float fluxFirst = 0;
+  float fluxLast = 0;
+  for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) {
+    if (i <= firstBelow) {
+      fluxFirst += flux->data.F32[i];
+    }
+    if (i <= lastAbove) {
+      fluxLast += flux->data.F32[i];
+    }
+  }
+  float fluxRPSum    = 0.5*(fluxLast + fluxFirst);
+  float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst);
+  // XXX need to use the weight appropriately here...
+
+  float rad     = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]);
+  float radErr  = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]);
+
+  if (!source->extpars->petrosian) {
+    source->extpars->petrosian = pmSourcePetrosianValuesAlloc ();
+  }
+  
+  // these are uncalibrated: instrumental mags and pixel units
+  source->extpars->petrosian->mag    = -2.5*log10(fluxRPSum);
+  source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum;
+
+  source->extpars->petrosian->rad    = rad;
+  source->extpars->petrosian->radErr = radErr;
+
   return true;
 

trunk/psphot/src/psphotRadialProfile.c

@@ -5 +5 @@
     // allocate pmSourceExtendedParameters, if not already defined
     if (!source->extpars) {
-        source->extpars = pmSourceExtendedParametersAlloc ();
+        source->extpars = pmSourceExtendedParsAlloc ();
     }
 
     if (!source->extpars->profile) {
-        source->extpars->profile = pmSourceRadialProfile (); 
+        source->extpars->profile = pmSourceRadialProfileAlloc (); 
     }    
     
@@ -19 +19 @@
     psVector *radius = source->extpars->profile->radius;
     psVector *flux   = source->extpars->profile->flux;
-    psVector *weight = source->extpars->profile->radius;
+    psVector *weight = source->extpars->profile->weight;
 
     // XXX use the extended source model here for Xo, Yo?
@@ -40 +40 @@
     flux->n = n;
 
-    SortVectorsByRadius (radius, flux, weight);
+    // XXX need to sort here
+    // SortVectorsByRadius (radius, flux, weight);
 
     return true;