IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5607


Ignore:
Timestamp:
Nov 25, 2005, 5:00:41 PM (20 years ago)
Author:
eugene
Message:

updating to work with current psModules

Location:
trunk/psphot
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r5593 r5607  
    1212IPSLIB  :=      $(shell pslib-config --cflags)
    1313
    14 INCS    =       $(IPSLIB)
    15 # LIBS  =       -lpsmodule $(LPSLIB)
    16 LIBS    =       $(LPSLIB)
     14LPSMOD  :=      $(shell psmodule-config --libs)
     15IPSMOD  :=      $(shell psmodule-config --cflags)
     16
     17INCS    =       $(IPSLIB) $(IPSMOD)
     18LIBS    =       $(LPSLIB) $(LPSMOD)
    1719CFLAGS  =       $(INCS) -std=c99 -Wall -Werror -g
    18 # CFLAGS        =       $(INCS) -std=c99 -g
    1920LFLAGS  =       $(LIBS)
    2021
    … …  
    3233$(SRC)/psphotSubtractPSF.$(ARCH).o \
    3334$(SRC)/psphotSortBySN.$(ARCH).o    \
    34 $(SRC)/psphotDefinePixels.$(ARCH).o \
     35$(SRC)/psphotDefinePixels.$(ARCH).o\
    3536$(SRC)/psphotMagnitudes.$(ARCH).o  \
    36 $(SRC)/psLibUtils.$(ARCH).o        \
    3737$(SRC)/psLine.$(ARCH).o            \
    38 $(SRC)/psEllipse.$(ARCH).o         \
    3938$(SRC)/psModulesUtils.$(ARCH).o    \
    4039$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
     40$(SRC)/psImageData.$(ARCH).o
     41
     42PSMODULES = \
     43$(SRC)/psEllipse.$(ARCH).o         \
    4144$(SRC)/pmPSF.$(ARCH).o             \
    4245$(SRC)/pmPSFtry.$(ARCH).o          \
    43 $(SRC)/psImageData.$(ARCH).o       \
    44 $(SRC)/pmModelGroup.$(ARCH).o       \
     46$(SRC)/pmModelGroup.$(ARCH).o      \
    4547$(SRC)/pmObjects_EAM.$(ARCH).o
    4648
    … …  
    6567
    6668psphot: $(BIN)/psphot.$(ARCH)
     69# $(BIN)/psphot.$(ARCH) : $(PSPHOT) $(PSMODULES)
    6770$(BIN)/psphot.$(ARCH) : $(PSPHOT)
    68 $(PSPHOT): $(SRC)/psphot.h
     71$(PSPHOT) $(PSMODULES): $(SRC)/psphot.h
    6972
    7073fitsource: $(BIN)/fitsource.$(ARCH)

  • trunk/psphot/doc/notes.txt

    r5593 r5607  
    11
    22Notes on psphot
     3
     42005.11.25
     5
     6  I've updated psphot to work with the current psLib v8, though a
     7  number of psLib fixes were needed.  These are pushed under
     8  psLib:eam_rel8_b2.  The psphot code which works with that version of
     9  psLib is tagged psphot_dev_04.
     10
     11  I'm working on converting psphot to work with the current release of
     12  psModules, which basically includes all of the object functions.
     13
     14  - change comments are relative to psModules (vs psphot versions)
     15
     16  psEllipse.c : OK
     17  psEllipse.h : OK
     18  pmModelGroup.c : OK (missing TGAUSS entry)
     19  pmModelGroup.h : OK (adds many comments)
     20
     21  pmPSFtry.c : a bunch of formatting weirdness (try-> got line breaks
     22               added pmPSFtryMetric_Alt
     23               fixed up usage of psVectorClipFitPolynomial1D 
     24               fixed up usage of psPolynomial..Alloc 
     25
     26  pmPSFtry.h : OK (Added extensive comments)
     27
     28  pmPSF.c : fixed up usage of psPolynomial..Alloc
     29
     30  pmPSF.h : OK (Added extensive comments)
     31
     32  pmObjects.c : cleaned up some formatting,
     33                fixed usage of psImageCountPixelMask
     34
     35  pmObjects.h : added apMag and fitMag to pmSource
     36                dropped PS_MODEL_name defines
     37
     38  psLibUtils.[ch] should be dropped from psModules
     39  psModulesUtils.[ch] should be dropped from psModules
     40
     41  Makefile.am : dropped psModulesUtils.[ch], psLibUtils.[ch]
     42
     43  - I have successfully tested psphot with the psModules.v8 code
     44    (eam_rel8_b1 tag).  I only needed to make a few modes in the
     45    psModules code. 
    346
    4472005.11.16

  • trunk/psphot/doc/versions.txt

    r5598 r5607  
    11
    222005.11.25
     3
     4  psphot_dev_05 should be compiled against psLib tag eam_rel8_b2 and
     5  psModules tag eam_rel8_b1.  this version removes all _EAM versions
     6  of psLib and psModules code.  It now uses the psModules version of
     7  the object code.
    38
    49  psphot_dev_04 should be compiled against psLib tag eam_rel8_b2 this

  • trunk/psphot/src/pmObjects_EAM.c

    r5593 r5607  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-11-25 23:13:43 $
     8 *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-11-26 03:00:41 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    723723    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    724724        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    725             if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
     725            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) {
     726                continue;
     727            }
    726728
    727729            psF32 xDiff = col + source->pixels->col0 - xPeak;
    … …  
    730732            // XXX EAM : calculate xDiff, yDiff up front;
    731733            //           radius is just a function of (xDiff, yDiff)
    732             if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
     734            if (!VALID_RADIUS(xDiff, yDiff, R2)) {
     735                continue;
     736            }
    733737
    734738            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
    … …  
    736740            // XXX EAM : check for valid S/N in pixel
    737741            // XXX EAM : should this limit be user-defined?
    738             if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) continue;
     742            if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) {
     743                continue;
     744            }
    739745           
    740746            Sum += pDiff;
    … …  
    10311037
    10321038        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1033         // allArray = psRegionSet (0, 0, 0, 0);
    10341039        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, allArray, PSPHOT_MASK_SATURATED);
    10351040
    … …  
    16701675        }
    16711676    }
    1672     # if (0)
    1673     else {
    1674         // try keeping the moments sky -- unsuccessful
    1675         for (int i = 1; i < paramMask->n; i++) {
    1676             paramMask->data.U8[i] = 0;
    1677         }
    1678         paramMask->data.U8[0] = 1;
    1679     }
    1680     # endif
    16811677
    16821678
    … …  
    18141810    return(p_pmSourceAddOrSubModel(image, mask, model, center, 1));
    18151811}
    1816 
    1817 
    1818 // XXX: Put this is psConstants.h
    1819 #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \
    1820 if (VEC1->n != N) { \
    1821     psError(PS_ERR_BAD_PARAMETER_SIZE, true, \
    1822             "psVector %s has size %d, should be %d.", \
    1823             #VEC1, VEC1->n, N); \
    1824     return(RVAL); \
    1825 }

  • trunk/psphot/src/pmObjects_EAM.h

    r5085 r5607  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-09-21 06:52:23 $
     7 *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-11-26 03:00:41 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2222#include<math.h>
    2323#include "pslib.h"
    24 #include "psLibUtils.h"
    25 // #include "pmModelGroup.h"
     24#include "psEllipse.h"
     25
     26# define psMemCopy(A)(psMemIncrRefCounter((A)))
    2627
    2728// defined mask values for mask image pixels.

  • trunk/psphot/src/pmPSFtry.c

    r5593 r5607  
    1616static void pmPSFtryFree (pmPSFtry *test) {
    1717
    18   if (test == NULL) return;
    19 
    20   psFree (test->psf);
    21   psFree (test->sources);
    22   psFree (test->modelFLT);
    23   psFree (test->modelPSF);
    24   psFree (test->metric);
    25   psFree (test->fitMag);
    26   psFree (test->mask);
    27   return;
     18    if (test == NULL) return;
     19
     20    psFree (test->psf);
     21    psFree (test->sources);
     22    psFree (test->modelFLT);
     23    psFree (test->modelPSF);
     24    psFree (test->metric);
     25    psFree (test->fitMag);
     26    psFree (test->mask);
     27    return;
    2828}
    2929
    … …  
    9595        // exclude the poor fits
    9696        if (!status) {
    97           try->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
    98           psFree (model);
    99           continue;
     97            try->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
     98            psFree (model);
     99            continue;
    100100        }
    101101        try->modelFLT->data[i] = model;
    … …  
    136136
    137137        // otherwise, save the resulting model
    138         // pmModelSkyOffset (modelPSF, x, y, RADIUS);
    139138        try->modelPSF->data[i] = modelPSF;
    140139
    … …  
    171170bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
    172171
    173   float dBin;
    174   int   nKeep, nSkip;
    175 
    176   // the measured (aperture - fit) magnitudes (dA == try->metric)
    177   //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    178   //     dA = dAo + dsky/flux
    179   //   where flux is the flux of the star
    180   // we fit this trend to find the infinite flux aperture correction (dAo),
    181   //   the nominal sky bias (dsky), and the error on dAo
    182   // the values of dA are contaminated by stars with close neighbors in the aperture
    183   //   we use an outlier rejection to avoid this bias
    184 
    185   FILE *f;
    186   f = fopen ("apresid.dat", "w");
    187   if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
    188 
    189   // rflux = ten(0.4*fitMag);
    190   psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    191   for (int i = 0; i < try->sources->n; i++) {
    192     if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    193     rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
    194     fprintf (f, "%3d %8.4f %12.5e %8.4f\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i]);
    195   }
    196   fclose (f);
    197 
    198   // XXX EAM : try 3hi/1lo sigma clipping on the rflux v dap fit
    199 
    200   // find min and max of (1/flux):
    201   psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    202   psVectorStats (stats, rflux, NULL, try->mask, PSFTRY_MASK_ALL);
     172    float dBin;
     173    int   nKeep, nSkip;
     174
     175    // the measured (aperture - fit) magnitudes (dA == try->metric)
     176    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     177    //     dA = dAo + dsky/flux
     178    //   where flux is the flux of the star
     179    // we fit this trend to find the infinite flux aperture correction (dAo),
     180    //   the nominal sky bias (dsky), and the error on dAo
     181    // the values of dA are contaminated by stars with close neighbors in the aperture
     182    //   we use an outlier rejection to avoid this bias
     183
     184    // rflux = ten(0.4*fitMag);
     185    psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     186    for (int i = 0; i < try->sources->n; i++) {
     187        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     188        rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
     189    }
     190
     191    // find min and max of (1/flux):
     192    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
     193    psVectorStats (stats, rflux, NULL, try->mask, PSFTRY_MASK_ALL);
    203194 
    204   // build binned versions of rflux, metric
    205   dBin = (stats->max - stats->min) / 10.0;
    206   psVector *rfBin = psVectorCreate (NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
    207   psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
    208   psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
    209   psFree (stats);
    210 
    211   psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
    212 
    213   // group data in daBin bins, measure lower 50% mean
    214   for (int i = 0; i < daBin->n; i++) {
    215 
    216     psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    217     tmp->n = 0;
    218 
    219     // accumulate data within bin range
    220     for (int j = 0; j < try->sources->n; j++) {
    221       // masked for: bad model fit, outlier in parameters
    222       if (try->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
     195    // build binned versions of rflux, metric
     196    dBin = (stats->max - stats->min) / 10.0;
     197    psVector *rfBin = psVectorCreate (NULL, stats->min, stats->max, dBin, PS_TYPE_F64);
     198    psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
     199    psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
     200    psFree (stats);
     201
     202    psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
     203
     204    // group data in daBin bins, measure lower 50% mean
     205    for (int i = 0; i < daBin->n; i++) {
     206
     207        psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     208        tmp->n = 0;
     209
     210        // accumulate data within bin range
     211        for (int j = 0; j < try->sources->n; j++) {
     212            // masked for: bad model fit, outlier in parameters
     213            if (try->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
    223214   
    224       // skip points with extreme dA values
    225       if (fabs(try->metric->data.F64[j]) > 0.5) continue;
    226 
    227       // skip points outside of this bin
    228       if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
    229       if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
    230 
    231       tmp->data.F64[tmp->n] = try->metric->data.F64[j];
    232       tmp->n ++;
    233     }
    234 
    235     // is this a valid point?
    236     maskB->data.U8[i] = 0;
    237     if (tmp->n < 2) {
    238       maskB->data.U8[i] = 1;
    239       psFree (tmp);
    240       continue;
    241     }
    242 
    243     // dA values are contaminated with low outliers
    244     // measure statistics only on upper 50% of points
    245     // this would be easier if we could sort in reverse:
    246 
    247     psVectorSort (tmp, tmp);
    248     nKeep = 0.5*tmp->n;
    249     nSkip = tmp->n - nKeep;
    250 
    251     psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
    252     for (int j = 0; j < tmp2->n; j++) {
    253       tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
    254     }
    255 
    256     stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    257     psVectorStats (stats, tmp2, NULL, NULL, 0);
    258     psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
    259 
    260     daBin->data.F64[i] = stats->sampleMedian;
    261 
     215            // skip points with extreme dA values
     216            if (fabs(try->metric->data.F64[j]) > 0.5) continue;
     217
     218            // skip points outside of this bin
     219            if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
     220            if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
     221
     222            tmp->data.F64[tmp->n] = try->metric->data.F64[j];
     223            tmp->n ++;
     224        }
     225
     226        // is this a valid point?
     227        maskB->data.U8[i] = 0;
     228        if (tmp->n < 2) {
     229            maskB->data.U8[i] = 1;
     230            psFree (tmp);
     231            continue;
     232        }
     233
     234        // dA values are contaminated with low outliers
     235        // measure statistics only on upper 50% of points
     236        // this would be easier if we could sort in reverse:
     237
     238        psVectorSort (tmp, tmp);
     239        nKeep = 0.5*tmp->n;
     240        nSkip = tmp->n - nKeep;
     241
     242        psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
     243        for (int j = 0; j < tmp2->n; j++) {
     244            tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
     245        }
     246
     247        stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
     248        psVectorStats (stats, tmp2, NULL, NULL, 0);
     249        psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
     250
     251        daBin->data.F64[i] = stats->sampleMedian;
     252
     253        psFree (stats);
     254        psFree (tmp);
     255        psFree (tmp2);
     256    }
     257
     258    // linear clipped fit to rfBin, daBin
     259    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     260    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     261    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
     262
     263    psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
     264    psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
     265
     266    stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
     267    stats = psVectorStats (stats, daResid, NULL, maskB, 1);
     268
     269    try->psf->ApResid = poly->coeff[0];
     270    try->psf->dApResid = stats->clippedStdev;
     271    try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     272
     273    psFree (rflux);
     274    psFree (rfBin);
     275    psFree (daBin);
     276    psFree (maskB);
     277    psFree (daBinFit);
     278    psFree (daResid);
     279    psFree (poly);
    262280    psFree (stats);
    263     psFree (tmp);
    264     psFree (tmp2);
    265   }
    266 
    267   // linear clipped fit to rfBin, daBin
    268   psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    269   psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    270   poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
    271 
    272   psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    273   psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    274 
    275   stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
    276   stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    277 
    278   try->psf->ApResid = poly->coeff[0];
    279   try->psf->dApResid = stats->clippedStdev;
    280   try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    281 
    282   psFree (rflux);
    283   psFree (rfBin);
    284   psFree (daBin);
    285   psFree (maskB);
    286   psFree (daBinFit);
    287   psFree (daResid);
    288   psFree (poly);
    289   psFree (stats);
    290   psFree (fitstat);
    291 
    292   return true;
     281    psFree (fitstat);
     282
     283    return true;
    293284}
    294285
    295286bool pmPSFtryMetric_Alt (pmPSFtry *try, float RADIUS) {
    296287
    297   // the measured (aperture - fit) magnitudes (dA == try->metric)
    298   //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    299   //     dA = dAo + dsky/flux
    300   //   where flux is the flux of the star
    301   // we fit this trend to find the infinite flux aperture correction (dAo),
    302   //   the nominal sky bias (dsky), and the error on dAo
    303   // the values of dA are contaminated by stars with close neighbors in the aperture
    304   //   we use an outlier rejection to avoid this bias
    305 
    306   // rflux = ten(0.4*fitMag);
    307   psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    308   for (int i = 0; i < try->sources->n; i++) {
    309     if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
    310     rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
    311   }
    312 
    313   // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
    314   psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    315 
    316   // XXX EAM
    317   stats->min = 1.0;
    318   stats->max = 3.0;
    319   stats->clipIter = 3;
    320 
    321   // linear clipped fit to rfBin, daBin
    322   psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
    323   poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
    324   fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
    325 
    326   try->psf->ApResid = poly->coeff[0];
    327   try->psf->dApResid = stats->sampleStdev;
    328   try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    329 
    330   FILE *f;
    331   f = fopen ("apresid.dat", "w");
    332   if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
    333 
    334   for (int i = 0; i < try->sources->n; i++) {
    335     fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
    336   }
    337   fclose (f);
    338 
    339   psFree (rflux);
    340   psFree (poly);
    341   psFree (stats);
    342 
    343   // psFree (daFit);
    344   // psFree (daResid);
    345 
    346   return true;
    347 }
     288    // the measured (aperture - fit) magnitudes (dA == try->metric)
     289    //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     290    //     dA = dAo + dsky/flux
     291    //   where flux is the flux of the star
     292    // we fit this trend to find the infinite flux aperture correction (dAo),
     293    //   the nominal sky bias (dsky), and the error on dAo
     294    // the values of dA are contaminated by stars with close neighbors in the aperture
     295    //   we use an outlier rejection to avoid this bias
     296
     297    // rflux = ten(0.4*fitMag);
     298    psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     299    for (int i = 0; i < try->sources->n; i++) {
     300        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     301        rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
     302    }
     303
     304    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
     305    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
     306
     307    // XXX EAM
     308    stats->min = 1.0;
     309    stats->max = 3.0;
     310    stats->clipIter = 3;
     311
     312    // linear clipped fit to rfBin, daBin
     313    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     314    poly = psVectorClipFitPolynomial1D (poly, stats, try->mask, PSFTRY_MASK_ALL, try->metric, NULL, rflux);
     315    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
     316
     317    try->psf->ApResid = poly->coeff[0];
     318    try->psf->dApResid = stats->sampleStdev;
     319    try->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     320
     321    FILE *f;
     322    f = fopen ("apresid.dat", "w");
     323    if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
     324
     325    for (int i = 0; i < try->sources->n; i++) {
     326        fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i], try->mask->data.U8[i]);
     327    }
     328    fclose (f);
     329
     330    psFree (rflux);
     331    psFree (poly);
     332    psFree (stats);
     333
     334    // psFree (daFit);
     335    // psFree (daResid);
     336
     337    return true;
     338}

  • trunk/psphot/src/psLibUtils.c

    r5593 r5607  
    6464    return (TRUE);
    6565}
    66 
    67 // XXX EAM a utility function
    68 bool p_psVectorPrintRow (int fd, psVector *a, char *name)
    69 {
    70 
    71     char line[1024];
    72 
    73     sprintf (line, "vector: %s\n", name);
    74     write (fd, line, strlen(line));
    75 
    76     for (int i = 0; i < a[0].n; i++) {
    77         sprintf (line, "%f  ", p_psVectorGetElementF64(a, i));
    78         write (fd, line, strlen(line));
    79     }
    80     write (fd, "\n", 1);
    81     return (true);
    82 }

  • trunk/psphot/src/psLibUtils.h

    r5351 r5607  
    22# ifndef PS_LIB_UTILS
    33# define PS_LIB_UTILS
    4 
    5 // XXX EAM : psEllipse needs to be be included in psLib
    6 # include "psEllipse.h"
    74
    85// structure to carry a dynamic string
    … …  
    1310} psLine;
    1411
    15 # define psMemCopy(A)(psMemIncrRefCounter((A)))
    16 
    17 // XXX EAM : my version using varience instead of stdev
    18 bool psMinimizeGaussNewtonDelta_EAM (psVector *delta,
    19                                      const psVector *params,
    20                                      const psVector *paramMask,
    21                                      const psArray  *x,
    22                                      const psVector *y,
    23                                      const psVector *yErr,
    24                                      psMinimizeLMChi2Func func);
    25 
    26 // minimize : using varience vs sigma and parameter limits
    27 psBool       p_psMinLM_GuessABP_EAM (psImage  *Alpha, psVector *Beta, psVector *Params, const psImage  *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda);
    28 psBool       psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);
    29 psF64        p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda);
    30 
    31 int          psImageCountPixelMask_EAM (psImage *mask, psU8 value); // added to SDRS
    32 
    33 // psLib extra utilities
    34 # if (0) // added to psLib v8
    35 int          psArgumentGet (int argc, char **argv, char *arg); // added to SDRS
    36 int          psArgumentRemove (int N, int *argc, char **argv); // added to SDRS
    37 bool         psTimerClear (char *name);   // added to SDRS
    38 psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type); // added to SDRS
    39 # endif
    40 
    41 bool         psTimerStart (char *name);   // added to SDRS
    42 psF64        psTimerMark (char *name);    // added to SDRS
    43 void         psTimerFree ();              // added to SDRS (as psTimerStop)
    44 psS32        psLogArguments (int *argc, char **argv);   // added to SDRS (part of psArgumentVerbosity)
    45 psS32        psTraceArguments (int *argc, char **argv); // added to SDRS (part of psArgumentVerbosity)
    46 
    47 // basic image functions
    48 bool         psImageInit (psImage *image,...); // added to SDRS (v.15)
    49 void         psImageSmooth_EAM (psImage *image, float sigma, float Nsigma); // added to SDRS (v.15)
    50 
    5112// psLine functions -- keep out for now?
    5213psLine      *psLineAlloc (int Nline);
    … …  
    5415bool         psLineAdd (psLine *line, char *format, ...);
    5516
    56 // not included in the .h file -- these are fairly lame, keep out?
    57 bool p_psVectorPrint (int fd, psVector *a, char *name);
    58 bool p_psVectorPrintRow (int fd, psVector *a, char *name);
    59 
    60 // polynomial functions -- all in SDRS, not done for v.15
    61 // XXX the pslib 0.7.0 polynomials are *still* nTerm, not nOrder!!!
    62 psPolynomial1D *Polynomial1DAlloc_EAM(psPolynomialType type, psS32 nOrder);
    63 void            Polynomial1DDump_EAM(psPolynomial1D *poly);
    64 psF64           Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x);
    65 psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
    66 psPolynomial1D *VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly, psVector* mask, psMaskType maskValue, const psVector* y, const psVector* yErr, const psVector* x);
    67 psPolynomial1D *VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly, psStats *stats, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x);
    68 
    69 psPolynomial2D *Polynomial2DAlloc_EAM(psPolynomialType type, psS32 nXorder, psS32 nYorder);
    70 void            Polynomial2DDump_EAM(psPolynomial2D *poly);
    71 psF64           Polynomial2DEval_EAM(const psPolynomial2D* myPoly, psF64 x, psF64 y);
    72 psVector       *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
    73 psVector       *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
    74 psPolynomial2D *VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x, const psVector* y);
    75 psPolynomial2D *VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly, psStats *stats, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x, const psVector* y);
    76 
    7717# endif

  • trunk/psphot/src/psphotArguments.c

    r5528 r5607  
    1212  fprintf (stderr, "starting... %s\n", psLibVersion());
    1313  psLogSetFormat ("M");
    14   psLogArguments (argc, argv);
    15   psTraceArguments (argc, argv);
     14  psArgumentVerbosity (argc, argv);
    1615
    1716  // optional mask image - add to config
Note: See TracChangeset for help on using the changeset viewer.