IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 4977


Ignore:
Timestamp:
Sep 7, 2005, 6:32:40 PM (21 years ago)
Author:
eugene
Message:

cleanup, organization

Location:
trunk/psphot
Files:
9 added
1 deleted
22 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/Makefile

    r4954 r4977  
    3232$(SRC)/psphotSubtractPSF.$(ARCH).o \
    3333$(SRC)/psphotSortBySN.$(ARCH).o    \
    34 $(SRC)/pspsf.$(ARCH).o             \
     34$(SRC)/psphotDefinePixels.$(ARCH).o    \
    3535$(SRC)/psLibUtils.$(ARCH).o        \
     36$(SRC)/psLine.$(ARCH).o            \
     37$(SRC)/psMinimize.$(ARCH).o        \
     38$(SRC)/psPolynomials.$(ARCH).o     \
     39$(SRC)/psEllipse.$(ARCH).o         \
    3640$(SRC)/psModulesUtils.$(ARCH).o    \
    37 $(SRC)/pmSourceUtils.$(ARCH).o     \
    3841$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
    39 $(SRC)/psLine.$(ARCH).o            \
    40 $(SRC)/psPolynomials.$(ARCH).o     \
    41 $(SRC)/psMinimize.$(ARCH).o        \
     42$(SRC)/pmPSF.$(ARCH).o             \
     43$(SRC)/pmPSFtry.$(ARCH).o          \
    4244$(SRC)/psImageData.$(ARCH).o       \
    43 $(SRC)/psEllipse.$(ARCH).o         \
    44 $(SRC)/pmModelInit.$(ARCH).o       \
     45$(SRC)/pmModelGroup.$(ARCH).o       \
    4546$(SRC)/pmObjects_EAM.$(ARCH).o
     47
     48MODELS = \
     49$(SRC)/models/pmModel_GAUSS.c      \
     50$(SRC)/models/pmModel_PGAUSS.c     \
     51$(SRC)/models/pmModel_QGAUSS.c     \
     52$(SRC)/models/pmModel_SGAUSS.c     
     53
     54$(SRC)/pmModelGroup.$(ARCH).o: $(MODELS)
     55
     56# deprecated
    4657
    4758FITSOURCE = \

  • trunk/psphot/doc/notes.txt

    r4954 r4977  
    4747  - fixed 'center' option in pmSourceAdd/Sub
    4848
     49  lingering concerns:
     50  - Polynomials and related functions are STILL messed up.  this is a
     51    major priority!!!
    4952  - pmSourceFitModel calculates sqrt(var), psMinimize calculates sq(sqrt(var))
    50   - pmSourceContour is using the model flux, not the image flux
    5153  - does psImageAlloc zero the image?  if so, drop init routine
    5254    from pmSourcePSFClump
    … …  
    6062    the psImageMaskRegion, etc functions.
    6163  - pmSourceSetPixelsCircle uses a hard-coded mask value
    62   -
     64  - pmSourceContour is using the model flux, not the image flux
    6365
    64662005.09.05

  • trunk/psphot/doc/psphot.txt

    r4954 r4977  
     1
     2TBD priorities:
     3
     4- test psImageSmooth / remove psImageSmooth_EAM
     5- unify pmSourceLocalSky / pmSourceLocalSky_EAM
     6- unify pmSourceMoments_EAM / pmSourceMoments_EAM
     7- merge with image data hierarchy stuff
     8- better name for psPSF_Test
    19
    210Defined APIs:
    311
    4 pmMoments           *pmMomentsAlloc();
    5 pmModel             *pmModelAlloc();
    6 pmSource            *pmSourceAlloc();
    7 psVector            *pmFindVectorPeaks()
    8 psArray             *pmFindImagePeaks()
    9 psList              *pmCullPeaks()
    10 pmSource            *pmSourceLocalSky()
    11 bool                 pmSourceMoments()
    12 pmPSFClump           pmSourcePSFClump()
    13 bool                 pmSourceRoughClass()
    14 bool                 pmSourceSetPixelsCircle()
    15 pmModel             *pmSourceModelGuess()
    16 psArray             *pmSourceContour()
    17 bool                 pmSourceFitModel()
    18 bool                 pmSourceAddModel()
    19 bool                 pmSourceSubModel()
    20 int                  pmModelParameterCount ();
    21 char                *pmModelGetType ();
    22 pmModelType          pmModelSetType ();
    23 pmModelFunc          pmModelFunc_GetFunction ();
    24 pmModelFlux          pmModelFlux_GetFunction ();
    25 pmModelRadius        pmModelRadius_GetFunction ();
    26 pmModelLimits        pmModelLimits_GetFunction ();
    27 pmModelGuessFunc     pmModelGuessFunc_GetFunction ();
    28 pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction ();
    29 pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction ();
     12top-level psphot functions : these do not need to be in psModules
     13  psMetadata  *psphotArguments ();
     14  psImageData *psphotSetup ();
     15  psStats     *psphotImageStats ();
     16  psArray     *pmPeaksSigmaLimit ();
     17  psArray     *psphotSourceStats ();
     18  pmPSF       *psphotChoosePSF ();
     19  bool         psphotApplyPSF ();
     20  bool         psphotFitGalaxies ();
     21  void         psphotOutput ();
     22
     23psphot-specific functions
     24  bool         psphotMarkPSF ();
     25  bool         psphotSubtractPSF ();
     26  int          psphotSortBySN ();
     27  int          psphotSaveImage ();
     28
     29** add the following entries to the psLibSDRS
     30
     31psLib extra utilities
     32  bool         psTimerStart ();         - already in SDRS
     33  void         psTimerFree ();          - already in SDRS
     34  bool         psTimerClear ();         - add to psLibSDRS
     35  psF64        psTimerMark ();          - already in SDRS
     36  psS32        psLogArguments ();       - now: psArgumentVerbosity in SDRS
     37  psS32        psTraceArguments ();     - now: psArgumentVerbosity in SDRS
     38  int          psArgumentGet ();        - already in SDRS
     39  int          psArgumentRemove ();     - already in SDRS
     40  psVector    *psVectorCreate ();       - add to psLibSDRS
     41  int          psImageCountPixelMask ();- add to psLibSDRS
     42  psVector    *psGetRowVectorFromImage(); - add to psLibSDRS?
     43  bool         p_psVectorPrintRow ();     - add to psLibSDRS?
     44
     45// basic image functions
     46  bool         psImageInit ();          - already in SDRS
     47  void         psImageSmooth_EAM ();    - test psLib version and drop
     48
     49// psLine functions
     50  psLine      *psLineAlloc ();          - add to psLibSDRS
     51  bool         psLineInit ();           - add to psLibSDRS
     52  bool         psLineAdd ();            - add to psLibSDRS
     53
     54** compare with current implementation / update psLib
     55   we need to define a mechanism to carry the beta values and the
     56   parameter limits
     57
     58// minimize
     59  psBool       p_psMinLM_GuessABP_EAM ();
     60  psBool       psMinimizeLMChi2_EAM();
     61  psF64        p_psMinLM_dLinear ();  -- not included in pslib.h files
     62
     63** already in psLib : some need to be cleaned up by MHPCC
     64
     65// polynomial functions
     66  psF64           Polynomial2DEval();
     67  psImage        *psBuildSums2D();
     68  psPolynomial2D *VectorFitPolynomial2DOrd_EAM();
     69  psPolynomial2D *Polynomial2DAlloc();
     70  void            psPolynomial2DDump ();
     71  psPolynomial2D *RobustFit2D_nomask();
     72  psPolynomial2D *RobustFit2D();
     73  psVector       *Polynomial2DEvalVector();
     74
     75  psVector       *psBuildSums1D();
     76  void            psPolynomial1DDump ();
     77  psF64           Polynomial1DEval_EAM();
     78  psVector       *Polynomial1DEvalVector_EAM();
     79  psPolynomial1D *Polynomial1DAlloc();
     80  psPolynomial1D *VectorFitPolynomial1DOrd_EAM();
     81
     82** make sure these are all in the psModules SDRS:
     83
     84pmObjects.h:
     85  pmMoments           *pmMomentsAlloc();
     86  pmModel             *pmModelAlloc();
     87  pmSource            *pmSourceAlloc();
     88  psVector            *pmFindVectorPeaks()
     89  psArray             *pmFindImagePeaks()
     90  psList              *pmCullPeaks()
     91  pmSource            *pmSourceLocalSky()
     92  bool                 pmSourceMoments()
     93  pmPSFClump           pmSourcePSFClump()
     94  bool                 pmSourceRoughClass()
     95  bool                 pmSourceSetPixelsCircle()
     96  pmModel             *pmSourceModelGuess()
     97  psArray             *pmSourceContour()
     98  bool                 pmSourceFitModel()
     99  bool                 pmSourceAddModel()
     100  bool                 pmSourceSubModel()
     101  int                  pmModelParameterCount ();
     102  char                *pmModelGetType ();
     103  pmModelType          pmModelSetType ();
     104  pmModelFunc          pmModelFunc_GetFunction ();
     105  pmModelFlux          pmModelFlux_GetFunction ();
     106  pmModelRadius        pmModelRadius_GetFunction ();
     107  pmModelLimits        pmModelLimits_GetFunction ();
     108  pmModelGuessFunc     pmModelGuessFunc_GetFunction ();
     109  pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction ();
     110  pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction ();
    30111
    31112model function abstractions:
    … …  
    37118  typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
    38119  typedef bool (*pmModelFitStatusFunc)(pmModel *model);
     120
     121conversions between elliptical shape representations
     122  EllipseAxes  EllipseMomentsToAxes ();
     123  EllipseShape EllipseAxesToShape ();
     124  EllipseAxes  EllipseShapeToAxes ();
     125
     126output functions
     127  bool         pmSourcesWriteText ();   -- be careful about psImageData...
     128  bool         pmSourcesWriteOBJ ();    -- be careful about psImageData...
     129  bool         pmSourcesWriteCMP ();    -- be careful about psImageData...
     130  bool         pmSourcesWriteCMF ();    -- be careful about psImageData...
     131  bool         pmSourcesWriteSX ();     -- be careful about psImageData...
     132  int          pmSourcesDophotType ();
     133  bool         pmPeaksWriteText ();
     134  bool         pmMomentsWriteText ();
     135  bool         pmModelWritePSFs ();
     136  bool         pmModelWriteFLTs ();
     137  bool         pmModelWriteNULLs ();
     138
     139// psf utilities
     140  pmPSF       *pmPSFAlloc ();
     141  pmPSF_Test  *pmPSF_TestAlloc ();
     142  pmPSF_Test  *pmPSF_TestModel ();
     143  bool         pmPSFFromModels ();
     144  pmModel     *pmModelFromPSF ();
     145  bool         pmSourcePhotometry ();
     146  bool         pmPSFMetricModel ();
     147
     148// psModule extra utilities
     149  bool         pmSourceFitModel_EAM();  -- requires updated version of psMinimize...
     150  bool         pmSourceLocalSky_EAM (); -- compare with pmSourceLocalSky & choose
     151  bool         pmSourceMoments_EAM();   -- compare with pmSourceMoments & choose
     152  bool         pmSourceDefinePixels();  -- be careful about psImageData...
     153  bool         pmModelFitStatus ();
     154  int          pmSourceDophotType ();
     155  psF32        pmConfigLookupF32 ();    -- merge with Paul's image container hierarchy
     156
     157** this needs to be confronted with the phase2 / image container hierarchy
     158
     159// psImageData functions
     160  psImageData *psImageDataAlloc ();

  • trunk/psphot/src/models/pmModel_GAUSS.c

    r4954 r4977  
    124124    for (int i = 4; i < 7; i++) {
    125125        psPolynomial2D *poly = psf->params->data[i-4];
    126         out[i] = Polynomial2DEval (poly, out[2], out[3]);
     126        out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
    127127    }
    128128    return(true);

  • trunk/psphot/src/models/pmModel_PGAUSS.c

    r4954 r4977  
    140140    for (int i = 4; i < 7; i++) {
    141141        psPolynomial2D *poly = psf->params->data[i-4];
    142         out[i] = Polynomial2DEval (poly, out[2], out[3]);
     142        out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
    143143    }
    144144    return(true);

  • trunk/psphot/src/models/pmModel_QGAUSS.c

    r4954 r4977  
    165165    for (int i = 4; i < 8; i++) {
    166166        psPolynomial2D *poly = psf->params->data[i-4];
    167         out[i] = Polynomial2DEval (poly, out[2], out[3]);
     167        out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
    168168    }
    169169    return(true);

  • trunk/psphot/src/models/pmModel_SGAUSS.c

    r4954 r4977  
    300300    for (int i = 4; i < 9; i++) {
    301301        psPolynomial2D *poly = psf->params->data[i-4];
    302         out[i] = Polynomial2DEval (poly, out[2], out[3]);
     302        out[i] = Polynomial2DEval_EAM (poly, out[2], out[3]);
    303303    }
    304304    return(true);

  • trunk/psphot/src/pmObjects_EAM.c

    r4954 r4977  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2005-09-07 03:33:01 $
     8 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2005-09-08 04:32:40 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    1919#include "psConstants.h"
    2020#include "pmObjects_EAM.h"
     21#include "pmModelGroup.h"
     22
    2123
    2224/******************************************************************************
    … …  
    653655
    654656    //
    655     // Loop through the subimage mask, initialize mask to 0x01 (invalid pixel)
    656     // XXX EAM : use PSPHOT_MASK_INVALID?
     657    // Loop through the subimage mask & initialize PSPHOT_MASK_INVALID
    657658    //
    658659    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
    659660        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
    660             subImageMask->data.U8[row][col] = 0x01;
     661            subImageMask->data.U8[row][col] = PSPHOT_MASK_INVALID;
    661662        }
    662663    }
    … …  
    674675    for (psS32 row = StartRow; row < EndRow; row++) {
    675676        for (psS32 col = StartCol; col < EndCol; col++) {
    676             subImageMask->data.U8[row][col] = 0;
     677            subImageMask->data.U8[row][col] = PSPHOT_MASK_CLEAR;
    677678        }
    678679    }
    … …  
    11161117
    11171118        // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1118 # define MASK_SATURATED 0x02
    1119         int Nsatpix = psImageCountPixelMask (tmpSrc->mask, MASK_SATURATED);
     1119        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, PSPHOT_MASK_SATURATED);
    11201120
    11211121        // saturated star (size consistent with PSF or larger)
    … …  
    16151615    }       
    16161616
    1617     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     1617    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    16181618
    16191619    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    … …  
    16211621    for (int i = 0; i < dparams->n; i++) {
    16221622        if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
    1623         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     1623        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    16241624    }
    16251625 
    … …  
    16401640    // XXX EAM get the Gauss-Newton distance for fixed model parameters
    16411641    if (paramMask != NULL) {
    1642         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     1642        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
    16431643        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
    16441644        for (int i = 0; i < dparams->n; i++) {
    16451645            if (!paramMask->data.U8[i]) continue;
    1646             dparams->data.F32[i] = delta->data.F64[i];
     1646            dparams->data.F32[i] = delta->data.F32[i];
    16471647        }
    16481648        psFree (delta);

  • trunk/psphot/src/pmObjects_EAM.h

    r4954 r4977  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-09-07 03:33:01 $
     7 *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-09-08 04:32:40 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2222#include<math.h>
    2323#include "pslib.h"
    24 
    25 // XXX EAM : this function should be added to psLib
    26 int          psImageCountPixelMask (psImage *mask, psU8 value);
    27 
    28 // XXX EAM : bug : missing from psLib/*.h
    29 bool psMinimizeGaussNewtonDelta (psVector *delta,
    30                                  const psVector *params,
    31                                  const psVector *paramMask,
    32                                  const psArray  *x,
    33                                  const psVector *y,
    34                                  const psVector *yErr,
    35                                  psMinimizeLMChi2Func func);
    36 
    37 // XXX EAM : switch to psPoly when ready
    38 psF64           Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
    39 
    40 // XXX EAM : psEllipse needs to be be included in psLib
    41 # include "psEllipse.h"
     24#include "psLibUtils.h"
     25// #include "pmModelGroup.h"
     26
     27// defined mask values for mask image pixels.
     28// this is probably a bad solution: we will want to set mask values
     29// outside of the PSPHOT code.  Perhaps we can set up a registered
     30// set of mask values with specific meanings that other functions
     31// can add to or define?
     32enum {
     33    PSPHOT_MASK_CLEAR     = 0x00,
     34    PSPHOT_MASK_INVALID   = 0x01,
     35    PSPHOT_MASK_SATURATED = 0x02,
     36    PSPHOT_MASK_MARKED    = 0x08,
     37} psphotMaskValues;
    4238
    4339/** pmPeakType
    … …  
    9187} pmMoments;
    9288
    93 typedef int pmModelType;
    94 
    9589/** pmPSFClump data structure
    9690 * 
    … …  
    10599    float dY;
    106100} pmPSFClump;
     101
     102typedef int pmModelType;
    107103
    108104/** pmModel data structure
    … …  
    169165pmSource;
    170166
    171 /** pmPSF data structure
    172  * 
    173  * 
    174  * 
    175  */
    176 typedef struct
    177 {
    178     pmModelType type;                   ///< PSF Model in use
    179     psArray *params;                    ///< Model parameters (psPolynomial2D)
    180     float chisq;                        ///< PSF goodness statistic
    181     int nPSFstars;                      ///< number of stars used to measure PSF
    182 }
    183 pmPSF;
    184 
    185167pmPeak *pmPeakAlloc(
    186168    int x,                              ///< Row-coordinate in image space
    … …  
    316298);
    317299
    318 /**
    319  *
    320  *  The object model functions are defined to allow for the flexible addition
    321  *  of new object models. Every object model, with parameters represented by
    322  *  pmModel, has an associated set of functions which provide necessary support
    323  *  operations. A set of abstract functions allow the programmer to select the
    324  *  approriate function or property for a specific named object model.
    325  *
    326  */
    327 
    328 /**
    329  *
    330  *  This function is the model chi-square minimization function for this model.
    331  *
    332  */
    333 typedef psMinimizeLMChi2Func pmModelFunc;
    334 
    335 
    336 /**
    337  *
    338  * This function returns the integrated flux for the given model parameters.
    339  */
    340 typedef psF64 (*pmModelFlux)(const psVector *params);
    341 
    342 
    343 /**
    344  *
    345  *  This function returns the radius at which the given model and parameters
    346  *  achieves the given flux.
    347  *
    348  */
    349 typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
    350 
    351 /**
    352  *
    353  *  This function sets the model parameter limits vectors for the given model
    354  *
    355  */
    356 typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
    357 
    358 /**
    359  *
    360  *  This function provides the model guess parameters based on the details of
    361  *   the given source.
    362  *
    363  */
    364 typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
    365 
    366 
    367 /**
    368  *
    369  *  This function constructs the PSF model for the given source based on the
    370  *  supplied psf and the FLT model for the object.
    371  *
    372  */
    373 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
    374 
    375 
    376 /**
    377  *
    378  *  This function returns the success / failure status of the given model fit
    379  *
    380  */
    381 typedef bool (*pmModelFitStatusFunc)(pmModel *model);
    382 
    383 /**
    384  *
    385  *  Each of the function types above has a corresponding function which returns
    386  *  the function given the model type:
    387  *
    388  */
    389 pmModelFunc          pmModelFunc_GetFunction (pmModelType type);
    390 pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
    391 pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
    392 pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
    393 pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
    394 pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
    395 pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
    396 
    397 // pmModelGroup utility functions
    398 int                  pmModelParameterCount (pmModelType type);
    399 char                *pmModelGetType (pmModelType type);
    400 pmModelType          pmModelSetType (char *name);
    401 
    402 // structure to carry model group functions
    403 typedef struct {
    404     char *name;
    405     int nParams;
    406     pmModelFunc          modelFunc;
    407     pmModelFlux          modelFlux;
    408     pmModelRadius        modelRadius;
    409     pmModelLimits        modelLimits;
    410     pmModelGuessFunc     modelGuessFunc;
    411     pmModelFromPSFFunc   modelFromPSFFunc;
    412     pmModelFitStatusFunc modelFitStatusFunc;
    413 } pmModelGroup;
    414 
    415300#endif

  • trunk/psphot/src/pmSourceUtils.c

    r4954 r4977  
    11# include "psphot.h"
    2 
    3 bool pmSourceDefinePixels(pmSource *mySource,
    4                           const psImageData *imdata,
    5                           psF32 x,
    6                           psF32 y,
    7                           psF32 Radius)
    8 {
    9     psRegion srcRegion;
    10 
    11     // Grab a subimage of the original image of size (2 * outerRadius).
    12     srcRegion = psRegionForSquare (x, y, Radius);
    13     srcRegion = psRegionForImage (imdata->image, &srcRegion);
    14 
    15     mySource->pixels = psImageSubset(imdata->image, srcRegion);
    16     mySource->noise  = psImageSubset(imdata->noise, srcRegion);
    17     mySource->mask   = psImageSubset(imdata->mask,  srcRegion);
    18 
    19     return(mySource);
    20 }
    21 
    22 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {
    23 
    24     float obsSum = 0;
    25     float fitSum = 0;
    26     float sky = model->params->data.F32[0];
    27 
    28     pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    29     fitSum = modelFluxFunc (model->params);
    30 
    31     for (int ix = 0; ix < image->numCols; ix++) {
    32         for (int iy = 0; iy < image->numRows; iy++) {
    33             if (mask->data.U8[iy][ix]) continue;
    34             obsSum += image->data.F32[iy][ix] - sky;
    35         }
    36     }
    37     if (obsSum <= 0) return false;
    38     if (fitSum <= 0) return false;
    39 
    40     *fitMag = -2.5*log10(fitSum);
    41     *obsMag = -2.5*log10(obsSum);
    42     return (true);
    43 }
    44 

  • trunk/psphot/src/psLibUtils.c

    r4954 r4977  
    11# include <strings.h>  // for strncasecmp
    22# include <pslib.h>
     3# include "psLibUtils.h"
    34
    45// XXX EAM : this is NOT included in the C99 headers ??
    … …  
    157158    memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32));
    158159    return(tmpVector);
    159 }
    160 
    161 // extract config informatin from config data or from image header, if specified
    162 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
    163 
    164     char *source;
    165     char *keyword;
    166     psF32 value;
    167     psMetadataItem *item;
    168 
    169     // look for the entry in the config collection
    170     item = psMetadataLookup (config, name);
    171     if (item == NULL) {
    172         psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
    173         value = psMetadataLookupF32 (status, header, name);
    174         if (*status == false) {
    175             psAbort ("pmConfigLookupF32", "no key %s in header", name);
    176         }
    177         *status = true;
    178         return (value);
    179     }   
    180 
    181     // I'm either expecting a string, with the name "HD:keyword"...
    182     if (item->type == PS_META_STR) {
    183         source = item->data.V;
    184         if (!strncasecmp (source, "HD:", 3)) {
    185             keyword = &source[3];
    186             value = psMetadataLookupF32 (status, header, keyword);
    187             if (*status == false) {
    188                 psAbort ("pmConfigLookupF32", "no key %s in config", name);
    189             }
    190             *status = true;
    191             // psFree (item);
    192             return (value);
    193         }       
    194     }
    195 
    196     //  ... or a value (F32?)
    197     if (item->type == PS_META_F32) {
    198         value = item->data.F32;
    199         // psFree (item);
    200         return (value);
    201     }
    202 
    203     *status = false;
    204     psError(PS_ERR_UNKNOWN, true, "invalid item");
    205     return (0);
    206160}
    207161

  • trunk/psphot/src/psLine.c

    r4954 r4977  
    1 # include "psphot.h"
     1# include <pslib.h>
     2# include "psLibUtils.h"
    23
    34static void psLineFree (psLine *line) {

  • trunk/psphot/src/psMinimize.c

    r4954 r4977  
    1 # include "psphot.h"
     1# include <pslib.h>
     2# include "psLibUtils.h"
    23
    34// XXX EAM : this file contains my alternate implementations of psMinimizeLMChi2

  • trunk/psphot/src/psModulesUtils.c

    r4954 r4977  
    1 # include "psphot.h"
     1# include "psModulesUtils.h"
     2
     3// extract config informatin from config data or from image header, if specified
     4psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
     5
     6    char *source;
     7    char *keyword;
     8    psF32 value;
     9    psMetadataItem *item;
     10
     11    // look for the entry in the config collection
     12    item = psMetadataLookup (config, name);
     13    if (item == NULL) {
     14        psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
     15        value = psMetadataLookupF32 (status, header, name);
     16        if (*status == false) {
     17            psAbort ("pmConfigLookupF32", "no key %s in header", name);
     18        }
     19        *status = true;
     20        return (value);
     21    }   
     22
     23    // I'm either expecting a string, with the name "HD:keyword"...
     24    if (item->type == PS_META_STR) {
     25        source = item->data.V;
     26        if (!strncasecmp (source, "HD:", 3)) {
     27            keyword = &source[3];
     28            value = psMetadataLookupF32 (status, header, keyword);
     29            if (*status == false) {
     30                psAbort ("pmConfigLookupF32", "no key %s in config", name);
     31            }
     32            *status = true;
     33            // psFree (item);
     34            return (value);
     35        }       
     36    }
     37
     38    //  ... or a value (F32?)
     39    if (item->type == PS_META_F32) {
     40        value = item->data.F32;
     41        // psFree (item);
     42        return (value);
     43    }
     44
     45    *status = false;
     46    psError(PS_ERR_UNKNOWN, true, "invalid item");
     47    return (0);
     48}
    249
    350// XXX EAM : these are my alternate implementations of psModule functions
    … …  
    321368    return (status);
    322369}
     370
     371bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {
     372
     373    float obsSum = 0;
     374    float fitSum = 0;
     375    float sky = model->params->data.F32[0];
     376
     377    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
     378    fitSum = modelFluxFunc (model->params);
     379
     380    for (int ix = 0; ix < image->numCols; ix++) {
     381        for (int iy = 0; iy < image->numRows; iy++) {
     382            if (mask->data.U8[iy][ix]) continue;
     383            obsSum += image->data.F32[iy][ix] - sky;
     384        }
     385    }
     386    if (obsSum <= 0) return false;
     387    if (fitSum <= 0) return false;
     388
     389    *fitMag = -2.5*log10(fitSum);
     390    *obsMag = -2.5*log10(obsSum);
     391    return (true);
     392}
     393

  • trunk/psphot/src/psPolynomials.c

    r4954 r4977  
    1 # include "psphot.h"
     1# include <pslib.h>
     2# include "psLibUtils.h"
    23
    34// XXX EAM : this file defines alternate versions of the polynomial fitting
    … …  
    910
    1011// write out the terms of the given 1D polynomial
    11 void psPolynomial1DDump (psPolynomial1D *poly) {
     12void Polynomial1DDump_EAM (psPolynomial1D *poly) {
    1213
    1314    for (int i = 0; i < poly->n + 1; i++) {
    … …  
    1617}   
    1718
    18 psF64 Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly)
    19 {
    20     psS32 loop_x = 0;
    21     psF64 polySum = 0.0;
    22     psF64 xSum = 1.0;
    23 
    24     for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
    25         if (myPoly->mask[loop_x] == 0) {
    26             polySum += xSum * myPoly->coeff[loop_x];
    27         }
    28         xSum *= x;
    29     }
    30 
    31     return(polySum);
    32 }
    33 
    34 psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
    35                                        const psVector *x)
    36 {
    37     // PS_POLY_CHECK_NULL(myPoly, NULL);
    38     // PS_VECTOR_CHECK_NULL(x, NULL);
    39     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    40 
    41     psVector *tmp;
    42 
    43     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    44     for (psS32 i=0;i<x->n;i++) {
    45         tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
    46     }
    47 
    48     return(tmp);
    49 }
    50 
    51 static void psPolynomial1DFree (psPolynomial1D *poly) {
    52 
    53   if (poly == NULL) return;
    54   psFree (poly->coeff);
    55   psFree (poly->coeffErr);
    56   psFree (poly->mask);
    57   return;
     19static void Polynomial1DFree_EAM (psPolynomial1D *poly) {
     20
     21    if (poly == NULL) return;
     22    psFree (poly->coeff);
     23    psFree (poly->coeffErr);
     24    psFree (poly->mask);
     25    return;
    5826}
    5927
    6028// XXX EAM : use Nterm = Norder + 1 definition 
    6129// XXX EAM : should we provide both order and nterms in struct?
    62 psPolynomial1D* Polynomial1DAlloc(psS32 nOrder,
    63                                   psPolynomialType type)
     30psPolynomial1D* Polynomial1DAlloc_EAM(psPolynomialType type,
     31                                      psS32 nOrder)
    6432{
    6533    // PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
    … …  
    8048        newPoly->mask[i] = 0;
    8149    }
    82     psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial1DFree);
     50    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial1DFree_EAM);
    8351    return(newPoly);
    8452}
    8553
     54psF64 Polynomial1DEval_EAM(const psPolynomial1D* myPoly, psF64 x)
     55{
     56    psS32 loop_x = 0;
     57    psF64 polySum = 0.0;
     58    psF64 xSum = 1.0;
     59
     60    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
     61        if (myPoly->mask[loop_x] == 0) {
     62            polySum += xSum * myPoly->coeff[loop_x];
     63        }
     64        xSum *= x;
     65    }
     66
     67    return(polySum);
     68}
     69
     70psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
     71                                     const psVector *x)
     72{
     73    // PS_POLY_CHECK_NULL(myPoly, NULL);
     74    // PS_VECTOR_CHECK_NULL(x, NULL);
     75    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     76
     77    psVector *tmp;
     78
     79    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
     80    for (psS32 i=0;i<x->n;i++) {
     81        tmp->data.F64[i] = Polynomial1DEval_EAM(myPoly, x->data.F64[i]);
     82    }
     83
     84    return(tmp);
     85}
     86
    8687// XXX EAM : my alternate BuildSums1D
    87 static psVector *BuildSums1D(psVector* sums,
    88                             psF64 x,
    89                             psS32 nTerm)
     88static psVector *BuildSums1D_EAM(psVector* sums,
     89                                psF64 x,
     90                                psS32 nTerm)
    9091{
    9192    psS32 nSum = 0;
    … …  
    109110
    110111// XXX EAM : test version of 1d fitting
    111 psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
    112                                              psVector *mask,
    113                                              const psVector *x,
    114                                              const psVector *y,
    115                                              const psVector *yErr)
     112psPolynomial1D* VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly,
     113                                          psVector *mask,
     114                                          psMaskType maskValue,
     115                                          const psVector *y,
     116                                          const psVector *yErr,
     117                                          const psVector *x)
    116118{
    117119    // I think this is 1 dimension down
    … …  
    130132    // XXX EAM : change from FILE to fd breaks this code:
    131133    if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
    132       FILE *f = psTraceGetDestination ();
    133       fprintf (f, "VectorFitPolynomial1D()\n");
    134       for (int i = 0; i < x->n; i++) {
    135         fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
    136       }
     134        FILE *f = psTraceGetDestination ();
     135        fprintf (f, "VectorFitPolynomial1D()\n");
     136        for (int i = 0; i < x->n; i++) {
     137            fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
     138        }
    137139    }
    138140# endif
    … …  
    155157    // Build the B and A data structs.
    156158    for (int k = 0; k < x->n; k++) {
    157         if ((mask != NULL) && mask->data.U8[k]) continue;
    158         xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
     159        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
     160        xSums = BuildSums1D_EAM(xSums, x->data.F64[k], nTerm);
    159161     
    160162        if (yErr == NULL) {
    … …  
    190192    }
    191193    else
    192     // LUD version of the fit
     194        // LUD version of the fit
    193195    {
    194196        psImage *ALUD = NULL;
    … …  
    218220// ********************** 2D polynomial functions ******************
    219221
    220 // XXX EAM : this version uses myPoly->nX as Norder, not Nterms
    221 psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
    222                                  const psVector *x,
    223                                  const psVector *y)
    224 
    225 {
    226     // PS_POLY_CHECK_NULL(myPoly, NULL);
    227     // PS_VECTOR_CHECK_NULL(x, NULL);
    228     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
    229     // PS_VECTOR_CHECK_NULL(y, NULL);
    230     // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
    231 
    232     psVector *tmp;
    233     psS32 vecLen=x->n;
    234 
    235     // Determine the length of the output vector to by the minimum of the x,y vectors
    236     if (y->n < vecLen) {
    237         vecLen = y->n;
    238     }
    239 
    240     // Create output vector to return
    241     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
    242 
    243     // Evaluate the polynomial at the specified points
    244     for (psS32 i=0; i<vecLen; i++) {
    245         tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
    246     }
    247 
    248     // Return output vector
    249     return(tmp);
    250 }
    251 
    252 // XXX EAM : this version uses the F64 vectors
    253 psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
    254                                   const psVector *x,
    255                                   const psVector *y)
    256 
    257 {
    258     // PS_POLY_CHECK_NULL(myPoly, NULL);
    259     // PS_VECTOR_CHECK_NULL(x, NULL);
    260     // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
    261     // PS_VECTOR_CHECK_NULL(y, NULL);
    262     // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
    263 
    264     psVector *tmp;
    265     psS32 vecLen=x->n;
    266 
    267     // Determine the length of the output vector to by the minimum of the x,y vectors
    268     if (y->n < vecLen) {
    269         vecLen = y->n;
    270     }
    271 
    272     // Create output vector to return
    273     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    274 
    275     // Evaluate the polynomial at the specified points
    276     for (psS32 i=0; i<vecLen; i++) {
    277         tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
    278     }
    279 
    280     // Return output vector
    281     return(tmp);
    282 }
    283 
    284 static void psPolynomial2DFree(psPolynomial2D *poly) {
    285 
    286   if (poly == NULL) return;
    287   for (int i = 0; i < poly->nX + 1; i++) {
    288     psFree (poly->coeff[i]);
    289     psFree (poly->coeffErr[i]);
    290     psFree (poly->mask[i]);
    291   }
    292   psFree (poly->coeff);
    293   psFree (poly->coeffErr);
    294   psFree (poly->mask);
    295   return;
     222static void Polynomial2DFree_EAM(psPolynomial2D *poly) {
     223
     224    if (poly == NULL) return;
     225    for (int i = 0; i < poly->nX + 1; i++) {
     226        psFree (poly->coeff[i]);
     227        psFree (poly->coeffErr[i]);
     228        psFree (poly->mask[i]);
     229    }
     230    psFree (poly->coeff);
     231    psFree (poly->coeffErr);
     232    psFree (poly->mask);
     233    return;
    296234}
    297235
    298236// XXX EAM : use Nterm = Norder + 1 definition 
    299237// the user requests a polynomial of order Norder
    300 psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
    301                                     psPolynomialType type)
     238psPolynomial2D* Polynomial2DAlloc_EAM(psPolynomialType type,
     239                                      psS32 nXorder,
     240                                      psS32 nYorder)
    302241{
    303242    // PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
    … …  
    332271        }
    333272    }
    334     psMemSetDeallocator(newPoly, (psFreeFunc) psPolynomial2DFree);
     273    psMemSetDeallocator(newPoly, (psFreeFunc) Polynomial2DFree_EAM);
    335274    return(newPoly);
    336275}
    337276
    338 // XXX EAM : BuildSums2D in analogy with BuildSums1D
    339 static psImage *BuildSums2D(psImage* sums,
    340                             psF64 x,      psF64 y,
    341                             psS32 nXterm, psS32 nYterm)
    342 {
    343     psS32 nXsum = 0;
    344     psS32 nYsum = 0;
    345     psF64 xSum = 1.0;
    346     psF64 ySum = 1.0;
    347 
    348     nXsum = 2*nXterm;
    349     nYsum = 2*nYterm;
    350     if (sums == NULL) {
    351         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    352     }
    353     if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
    354         psFree (sums);
    355         sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
    356     }
    357 
    358     ySum = 1.0;
    359     for (int j = 0; j < nYsum; j++) {
    360         xSum = ySum;
    361         for (int i = 0; i < nXsum; i++) {
    362             sums->data.F64[i][j] = xSum;
    363             xSum *= x;
    364         }
    365         ySum *= y;
    366     }
    367     return (sums);
    368 }
    369 
    370 // XXX EAM : test version of 2d fitting
    371 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
    372                                              psVector* mask,
    373                                              const psVector* x,
    374                                              const psVector* y,
    375                                              const psVector* z,
    376                                              const psVector* zErr)
    377 {
    378     // I think this is 1 dimension down
    379     psImage*      A = NULL;
    380     psVector*     B = NULL;
    381     psImage*   Sums = NULL;
    382     psF64 wt;
    383     psS32 nTerm, nXterm, nYterm;
    384 
    385     nXterm = myPoly->nX + 1;
    386     nYterm = myPoly->nY + 1;
    387     nTerm = nXterm * nYterm;
    388 
    389     A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
    390     B = psVectorAlloc(nTerm, PS_TYPE_F64);
    391 
    392     // Initialize data structures (why is this not a function!)
    393     for (int i = 0; i < nTerm; i++) {
    394         B->data.F64[i] = 0.0;
    395         for (int j = 0; j < nTerm; j++) {
    396             A->data.F64[i][j] = 0.0;
    397         }
    398     }
    399 
    400     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
    401 
    402     // Build the B and A data structs.
    403     for (int k  = 0; k < x->n; k++) {
    404       if ((mask != NULL) && mask->data.U8[k]) continue;
    405         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
    406      
    407         if (zErr == NULL) {
    408             wt = 1.0;
    409         } else {
    410             // this should probably by zErr^2 !!
    411             // this should filter zErr == 0 values
    412             wt = 1.0 / zErr->data.F64[k];
    413         }
    414 
    415         // we could skip half of the array and assign at the end
    416         // we must handle masked orders
    417         for (int n = 0; n < nXterm; n++) {
    418             for (int m = 0; m < nYterm; m++) {
    419                 B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
    420             }
    421         }
    422 
    423         for (int i = 0; i < nXterm; i++) {
    424             for (int j = 0; j < nYterm; j++) {
    425                 for (int n = 0; n < nXterm; n++) {
    426                     for (int m = 0; m < nYterm; m++) {
    427                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
    428                     }
    429                 }
    430             }
    431         }
    432     }
    433  
    434     // does the solution in place
    435     psGaussJordan (A, B);
    436    
    437     // XXX: How do we know if these routines were successful?
    438     // ALUD = psMatrixLUD(ALUD, &outPerm, A);
    439     // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
    440 
    441     for (int n = 0; n < nXterm; n++) {
    442         for (int m = 0; m < nYterm; m++) {
    443             myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
    444         }
    445     }
    446 
    447     psFree(A);
    448     psFree(B);
    449     psFree(Sums);
    450 
    451     psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
    452             "---- VectorFitPolynomial2DOrd() begin ----\n");
    453     return (myPoly);
    454 }
    455 
    456277// write out the terms of the given 2D polynomial
    457 void psPolynomial2DDump (psPolynomial2D *poly) {
     278void Polynomial2DDump_EAM (psPolynomial2D *poly) {
    458279
    459280    for (int i = 0; i < poly->nX + 1; i++) {
    … …  
    464285}   
    465286
    466 psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
    467                             const psVector* x,
    468                             const psVector* y,
    469                             const psVector* z,
    470                             const psVector* dz)
    471 {
    472     psVector *X;
    473     psVector *Y;
    474     psVector *Z;
    475     psVector *dZ;
    476 
    477     psVector *zFit   = NULL;
    478     psVector *zResid = NULL;
    479     psStats  *stats  = NULL;
    480 
    481     X  = psVectorCopy (NULL, x, PS_TYPE_F64);
    482     Y  = psVectorCopy (NULL, y, PS_TYPE_F64);
    483     Z  = psVectorCopy (NULL, z, PS_TYPE_F64);
    484     dZ = psVectorCopy (NULL, dz, PS_TYPE_F64);
    485 
    486     for (int N = 0; N < 3; N++) {
    487         // XXX EAM : this would be better defined with an element mask
    488         poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
    489         zFit   = Polynomial2DEvalVectorD (poly, x, y);
    490         zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
    491 
    492         stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    493         stats  = psVectorStats (stats, zResid, NULL, NULL, 0);
    494         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
    495 
    496         // re-create X, Y, Z, dZ if pts are valid
    497         int n = 0;
    498         for (int i = 0; i < zResid->n; i++) {
    499             if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
    500               continue;
    501             }
    502             X->data.F64[n]  =  x->data.F64[i];
    503             Y->data.F64[n]  =  y->data.F64[i];
    504             Z->data.F64[n]  =  z->data.F64[i];
    505             dZ->data.F64[n] = dz->data.F64[i];
    506             n++;
    507         }
    508         X->n = n;
    509         Y->n = n;
    510         Z->n = n;
    511         dZ->n = n;
    512     }
    513     return (poly);
    514 }
    515 
    516 // XXX EAM : be careful here with F32 vs F64 vectors
    517 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
    518                             psVector* mask,
    519                             const psVector* x,
    520                             const psVector* y,
    521                             const psVector* z,
    522                             const psVector* dz)
    523 {
    524     // PS_VECTOR_CHECK_NULL(mask, NULL);
    525     // PS_VECTOR_CHECK_NULL(x, NULL);
    526     // PS_VECTOR_CHECK_NULL(y, NULL);
    527     // PS_VECTOR_CHECK_NULL(z, NULL);
    528     // PS_VECTOR_CHECK_NULL(dz, NULL);
    529 
    530     psVector *zFit   = NULL;
    531     psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
    532     psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    533 
    534     for (int N = 0; N < 3; N++) {
    535         poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
    536         zFit   = Polynomial2DEvalVectorD (poly, x, y);
    537         zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
    538 
    539         stats  = psVectorStats (stats, zResid, NULL, mask, 1);
    540         psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
    541                  stats->sampleMean, stats->sampleStdev);
    542 
    543         // set mask if pts are not valid
    544         // we are masking out any point which is out of range
    545         // recovery is not allowed with this scheme
    546         for (int i = 0; i < zResid->n; i++) {
    547           if (mask->data.U8[i]) continue;
    548           if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
    549             mask->data.U8[i] = 1;
    550             continue;
    551           }       
    552         }
    553         psFree (zFit);
    554     }
    555     psFree (zResid);
    556     psFree (stats);
    557     return (poly);
    558 }
    559 
    560 // XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64)
    561 
    562 psF64 Polynomial2DEval(const psPolynomial2D* myPoly,
    563                        psF64 x,
    564                        psF64 y)
     287psF64 Polynomial2DEval_EAM(const psPolynomial2D* myPoly,
     288                           psF64 x,
     289                           psF64 y)
    565290{
    566291    // PS_POLY_CHECK_NULL(myPoly, NAN);
    … …  
    587312    return(polySum);
    588313}
     314
     315// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
     316psVector *Polynomial2DEvalVector_EAM(const psPolynomial2D *myPoly,
     317                                     const psVector *x,
     318                                     const psVector *y)
     319
     320{
     321    // PS_POLY_CHECK_NULL(myPoly, NULL);
     322    // PS_VECTOR_CHECK_NULL(x, NULL);
     323    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
     324    // PS_VECTOR_CHECK_NULL(y, NULL);
     325    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
     326
     327    psVector *tmp;
     328    psS32 vecLen=x->n;
     329
     330    // Determine the length of the output vector to by the minimum of the x,y vectors
     331    if (y->n < vecLen) {
     332        vecLen = y->n;
     333    }
     334
     335    // Create output vector to return
     336    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
     337
     338    // Evaluate the polynomial at the specified points
     339    for (psS32 i=0; i<vecLen; i++) {
     340        tmp->data.F32[i] = Polynomial2DEval_EAM(myPoly, x->data.F32[i], y->data.F32[i]);
     341    }
     342
     343    // Return output vector
     344    return(tmp);
     345}
     346
     347// XXX EAM : this version uses the F64 vectors
     348psVector *Polynomial2DEvalVectorD_EAM(const psPolynomial2D *myPoly,
     349                                  const psVector *x,
     350                                  const psVector *y)
     351
     352{
     353    // PS_POLY_CHECK_NULL(myPoly, NULL);
     354    // PS_VECTOR_CHECK_NULL(x, NULL);
     355    // PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
     356    // PS_VECTOR_CHECK_NULL(y, NULL);
     357    // PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
     358
     359    psVector *tmp;
     360    psS32 vecLen=x->n;
     361
     362    // Determine the length of the output vector to by the minimum of the x,y vectors
     363    if (y->n < vecLen) {
     364        vecLen = y->n;
     365    }
     366
     367    // Create output vector to return
     368    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
     369
     370    // Evaluate the polynomial at the specified points
     371    for (psS32 i=0; i<vecLen; i++) {
     372        tmp->data.F64[i] = Polynomial2DEval_EAM(myPoly,x->data.F64[i],y->data.F64[i]);
     373    }
     374
     375    // Return output vector
     376    return(tmp);
     377}
     378
     379// XXX EAM : BuildSums2D in analogy with BuildSums1D
     380static psImage *BuildSums2D_EAM(psImage* sums,
     381                                psF64 x,      psF64 y,
     382                                psS32 nXterm, psS32 nYterm)
     383{
     384    psS32 nXsum = 0;
     385    psS32 nYsum = 0;
     386    psF64 xSum = 1.0;
     387    psF64 ySum = 1.0;
     388
     389    nXsum = 2*nXterm;
     390    nYsum = 2*nYterm;
     391    if (sums == NULL) {
     392        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     393    }
     394    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
     395        psFree (sums);
     396        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
     397    }
     398
     399    ySum = 1.0;
     400    for (int j = 0; j < nYsum; j++) {
     401        xSum = ySum;
     402        for (int i = 0; i < nXsum; i++) {
     403            sums->data.F64[i][j] = xSum;
     404            xSum *= x;
     405        }
     406        ySum *= y;
     407    }
     408    return (sums);
     409}
     410
     411// XXX EAM : test version of 2d fitting
     412psPolynomial2D* VectorFitPolynomial2D_EAM(psPolynomial2D* myPoly,
     413                                             psVector* mask,
     414                                             psMaskType maskValue,
     415                                             const psVector* z,
     416                                             const psVector* zErr,
     417                                             const psVector* x,
     418                                             const psVector* y)
     419{
     420    // I think this is 1 dimension down
     421    psImage*      A = NULL;
     422    psVector*     B = NULL;
     423    psImage*   Sums = NULL;
     424    psF64 wt;
     425    psS32 nTerm, nXterm, nYterm;
     426
     427    nXterm = myPoly->nX + 1;
     428    nYterm = myPoly->nY + 1;
     429    nTerm = nXterm * nYterm;
     430
     431    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
     432    B = psVectorAlloc(nTerm, PS_TYPE_F64);
     433
     434    // Initialize data structures (why is this not a function!)
     435    for (int i = 0; i < nTerm; i++) {
     436        B->data.F64[i] = 0.0;
     437        for (int j = 0; j < nTerm; j++) {
     438            A->data.F64[i][j] = 0.0;
     439        }
     440    }
     441
     442    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
     443
     444    // Build the B and A data structs.
     445    for (int k  = 0; k < x->n; k++) {
     446        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) continue;
     447        Sums = BuildSums2D_EAM(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
     448     
     449        if (zErr == NULL) {
     450            wt = 1.0;
     451        } else {
     452            // this should probably by zErr^2 !!
     453            // this should filter zErr == 0 values
     454            wt = 1.0 / zErr->data.F64[k];
     455        }
     456
     457        // we could skip half of the array and assign at the end
     458        // we must handle masked orders
     459        for (int n = 0; n < nXterm; n++) {
     460            for (int m = 0; m < nYterm; m++) {
     461                B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
     462            }
     463        }
     464
     465        for (int i = 0; i < nXterm; i++) {
     466            for (int j = 0; j < nYterm; j++) {
     467                for (int n = 0; n < nXterm; n++) {
     468                    for (int m = 0; m < nYterm; m++) {
     469                        A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
     470                    }
     471                }
     472            }
     473        }
     474    }
     475 
     476    // does the solution in place
     477    psGaussJordan (A, B);
     478   
     479    // XXX: How do we know if these routines were successful?
     480    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
     481    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
     482
     483    for (int n = 0; n < nXterm; n++) {
     484        for (int m = 0; m < nYterm; m++) {
     485            myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
     486        }
     487    }
     488
     489    psFree(A);
     490    psFree(B);
     491    psFree(Sums);
     492
     493    psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
     494            "---- VectorFitPolynomial2DOrd() begin ----\n");
     495    return (myPoly);
     496}
     497
     498// XXX EAM : be careful here with F32 vs F64 vectors
     499psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
     500                                              psStats *stats,
     501                                              psVector* mask,
     502                                              psMaskType maskValue,
     503                                              const psVector* z,
     504                                              const psVector* dz,
     505                                              const psVector* x,
     506                                              const psVector* y)
     507{
     508    psVector *zFit   = NULL;
     509    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
     510
     511    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
     512    stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     513
     514    for (int N = 0; N < 3; N++) {
     515        poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
     516        zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
     517        zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
     518
     519        stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
     520        psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n",
     521                 stats->sampleMean, stats->sampleStdev);
     522
     523        // set mask if pts are not valid
     524        // we are masking out any point which is out of range
     525        // recovery is not allowed with this scheme
     526        for (int i = 0; i < zResid->n; i++) {
     527            if (mask->data.U8[i]) continue;
     528            if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
     529                mask->data.U8[i] &= 0x01;
     530                continue;
     531            }       
     532        }
     533        psFree (zFit);
     534    }
     535    psFree (zResid);
     536    return (poly);
     537}
     538
     539// XXX EAM : VectorFitPolynomial2D and Polynomial2DEvalVector require different types (F32 vs F64)

  • trunk/psphot/src/psphot.h

    r4954 r4977  
    33# include <unistd.h>   // for unlink
    44# include <pslib.h>
     5# include "psLibUtils.h"
    56# include "pmObjects_EAM.h"
     7# include "psModulesUtils.h"
     8# include "pmPSF.h"
     9# include "pmPSFtry.h"
     10# include "pmModelGroup.h"
    611
    712typedef struct {
    … …  
    1116    psMetadata *header;
    1217} psImageData;
    13 
    14 // data to test a given PSF model type
    15 typedef struct {
    16     pmModelType modelType;
    17     pmPSF      *psf;
    18     psArray    *sources;      // pointers to the original sources
    19     psArray    *modelFLT;     // model fits, floating parameters
    20     psArray    *modelPSF;     // model fits, PSF parameters
    21     psVector   *mask;
    22     psVector   *metric;
    23     psVector   *fitMag;
    24     float       ApResid;
    25     float       dApResid;
    26     float       skyBias;
    27 } pmPSF_Test;
    28 
    29 // structure to carry a dynamic string
    30 typedef struct {
    31     int NLINE;
    32     int Nline;
    33     char *line;
    34 } psLine;
    35 
    36 enum {
    37     PSPHOT_MASK_INVALID   = 0x01,
    38     PSPHOT_MASK_SATURATED = 0x02,
    39     PSPHOT_MASK_MARKED    = 0x08,
    40 } MaskValues;
    41 
    42 # define psMemCopy(A)(psMemIncrRefCounter((A)))
    4318
    4419// top-level psphot functions
    … …  
    5631bool         psphotSubtractPSF (pmSource *source);
    5732int          psphotSortBySN (const void **a, const void **b);
     33int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
     34bool         psphotDefinePixels (pmSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
    5835
    5936// psImageData functions
    6037psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header);
    61 
    62 // psf utilities
    63 pmPSF       *pmPSFAlloc (pmModelType type);
    64 pmPSF_Test  *pmPSF_TestAlloc (psArray *stars, char *modelName);
    65 pmPSF_Test  *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
    66 bool         pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
    67 pmModel     *pmModelFromPSF (pmModel *model, pmPSF *psf);
    68 bool         pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask);
    69 bool         pmPSFMetricModel (pmPSF_Test *test, float RADIUS);
    7038
    7139// output functions
    … …  
    8149bool         pmModelWriteFLTs (psArray *sources, char *filename);
    8250bool         pmModelWriteNULLs (psArray *sources, char *filename);
    83 
    84 // psModule extra utilities
    85 bool         pmSourceDefinePixels(pmSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
    86 bool         pmSourceLocalSky_EAM (pmSource *source, psStatsOptions statsOptions, psF32 Radius);
    87 bool         pmSourceFitModel_EAM(pmSource *source, pmModel *model, const bool PSF);
    88 bool         pmSourceMoments_EAM(pmSource *source, psF32 radius);
    89 bool         pmModelFitStatus (pmModel *model);
    90 int          pmSourceDophotType (pmSource *source);
    91 
    92 // minimize
    93 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);
    94 psBool       psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);
    95 psF64        p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda);
    96 
    97 // psLib extra utilities
    98 bool         psTimerStart (char *name);
    99 void         psTimerFree ();
    100 bool         psTimerClear (char *name);
    101 psF64        psTimerMark (char *name);
    102 psS32        psLogArguments (int *argc, char **argv);
    103 psS32        psTraceArguments (int *argc, char **argv);
    104 int          psArgumentGet (int argc, char **argv, char *arg);
    105 int          psArgumentRemove (int N, int *argc, char **argv);
    106 psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
    107 psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type);
    108 psVector    *psGetRowVectorFromImage(psImage *image, psU32 row);
    109 int          psImageCountPixelMask (psImage *mask, psU8 value);
    110 
    111 // basic image functions
    112 bool         psImageInit (psImage *image,...);
    113 void         psImageSmooth_EAM (psImage *image, float sigma, float Nsigma);
    114 int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
    115 
    116 // psLine functions
    117 psLine      *psLineAlloc (int Nline);
    118 bool         psLineInit (psLine *line);
    119 bool         psLineAdd (psLine *line, char *format, ...);
    120 
    121 // polynomial functions
    122 psF64           Polynomial2DEval(const psPolynomial2D* myPoly, psF64 x, psF64 y);
    123 psImage        *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
    124 psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    125 psPolynomial2D *Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
    126 void            psPolynomial2DDump (psPolynomial2D *poly);
    127 psPolynomial2D *RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    128 psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
    129 psVector       *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
    130 
    131 psVector       *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
    132 void            psPolynomial1DDump (psPolynomial1D *poly);
    133 psF64           Polynomial1DEval_EAM(psF64 x, const psPolynomial1D* myPoly);
    134 psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
    135 psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
    136 psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
    137 
    138 // XXX EAM : bug : missing from psLib/*.h
    139 bool psMinimizeGaussNewtonDelta (psVector *delta,
    140                                  const psVector *params,
    141                                  const psVector *paramMask,
    142                                  const psArray  *x,
    143                                  const psVector *y,
    144                                  const psVector *yErr,
    145                                  psMinimizeLMChi2Func func);
    146 
    147 // not included in the .h file?
    148 bool p_psVectorPrint (
    149     int fd,                            ///< output file descriptor
    150     psVector *a,                       ///< vector to print
    151     char *name                         ///< name of vector (for title)
    152 );
    153 
    154 bool p_psVectorPrintRow (int fd, psVector *a, char *name);

  • trunk/psphot/src/psphotApplyPSF.c

    r4954 r4977  
    6262        if (model->radius > OUTER_RADIUS) {
    6363          // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
    64           pmSourceDefinePixels (source, imdata, x, y, model->radius);
     64            psphotDefinePixels (source, imdata, x, y, model->radius);
    6565        }
    6666
     67        if (i > 66) psTraceSetLevel (".psLib.dataManip.psMinimizeLMChi2", 5);
     68           
    6769        // fit PSF model (set/unset the pixel mask)
    6870        psImageKeepCircle (source->mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);

  • trunk/psphot/src/psphotArguments.c

    r4954 r4977  
    99
    1010  // basic pslib options
    11   fprintf (stderr, "starting...\n");
     11  fprintf (stderr, "starting... %s\n", psLibVersion());
    1212  psLogSetFormat ("M");
    1313  psLogArguments (argc, argv);

  • trunk/psphot/src/psphotChoosePSF.c

    r4949 r4977  
    11# include "psphot.h"
    22
    3 // test PSF models and select best option
     3// try PSF models and select best option
    44
    55pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats)
    … …  
    77    bool            status;
    88    char           *modelName;
    9     pmPSF_Test     *test  = NULL;
     9    pmPSFtry       *try   = NULL;
    1010    psArray        *stars = NULL;
    1111    psMetadataItem *item  = NULL;
    … …  
    3434    psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE);
    3535
    36     // set up an array to store the test results
    37     psArray *tests = psArrayAlloc (list->n);
     36    // set up an array to store the results
     37    psArray *models = psArrayAlloc (list->n);
    3838
    39     // test each model option listed in config
    40     for (int i = 0; i < tests->n; i++) {
     39    // try each model option listed in config
     40    for (int i = 0; i < models->n; i++) {
    4141
    4242        item = psListGetAndIncrement (iter);
    4343        modelName = item->data.V;
    4444
    45         tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
     45        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
    4646        psFree (modelName);
    4747        psFree (item);
    … …  
    5151    psFree (stars);
    5252
    53     // select the best of the test models
     53    // select the best of the models
    5454    // here we are using the clippedStdev on the metric as the indicator
    55     test = tests->data[0];
     55    try = models->data[0];
    5656    int bestN = 0;
    57     float bestM = test->dApResid;
    58     for (int i = 1; i < tests->n; i++) {
    59         test = tests->data[i];
    60         float M = test->dApResid;
     57    float bestM = try->dApResid;
     58    for (int i = 1; i < models->n; i++) {
     59        try = models->data[i];
     60        float M = try->dApResid;
    6161        if (M < bestM) {
    6262            bestM = M;
    … …  
    6565    }
    6666
    67     // keep only the selected test:
    68     test = tests->data[bestN];
    69     modelName = pmModelGetType (test->modelType);
    70     psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid);
     67    // keep only the selected model:
     68    try = models->data[bestN];
     69    modelName = pmModelGetType (try->modelType);
     70    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, try->ApResid, try->dApResid);
    7171
    72     // XXX I am having trouble tracking down the double frees vs leaks
    73     # if (0)
    74     // we are keeping the modelPSF fits from the PSF test,
    75     // but we replace these later on (apply_psf_model)
    76     // set source->model based on best psf model fit
    77     for (int i = 0; i < test->sources->n; i++) {
    78         pmSource *source = test->sources->data[i];
    79         // drop masked sources from PSFSTAR list
    80         if (test->mask->data.U8[i]) {
    81           source->type = PS_SOURCE_OTHER;
    82           source->modelPSF = NULL;
    83         } else {
    84           // source->modelPSF = psMemCopy(test->modelPSF->data[i]);
    85           source->modelPSF = test->modelPSF->data[i];
    86         }
    87     }
    88     # endif
    89 
    90     pmPSF *psf = psMemCopy(test->psf);
    91     psFree (tests);
     72    // keep only the pmPSF resulting from this analysis
     73    pmPSF *psf = psMemCopy(try->psf);
     74    psFree (models);
    9275    return (psf);
    9376}

  • trunk/psphot/src/psphotFitGalaxies.c

    r4954 r4977  
    7272        if (model->radius > OUTER_RADIUS) {
    7373          // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
    74           pmSourceDefinePixels (source, imdata, x, y, model->radius);
     74          psphotDefinePixels (source, imdata, x, y, model->radius);
    7575        }
    7676

  • trunk/psphot/src/psphotSourceStats.c

    r4954 r4977  
    2424        // allocate image, noise, mask arrays for each peak (square of radius OUTER)
    2525        // XXX EAM : this has 12 leaks (v.5)
    26         pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
     26        psphotDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
    2727
    2828        // measure a local sky value

  • trunk/psphot/src/pspsf.c

    r4954 r4977  
    11# include "psphot.h"
    22
    3 static void pmPSFFree (pmPSF *psf) {
    4 
    5   if (psf == NULL) return;
    6 
    7   psFree (psf->params);
    8   return;
    9 }
    10 
    11 // a PSF always has 4 parameters fewer than the equivalent model
    12 pmPSF *pmPSFAlloc (pmModelType type) {
    13 
    14     int Nparams;
    15 
    16     pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
    17 
    18     psf->type = type;
    19     psf->chisq = 0.0;
    20 
    21     Nparams = pmModelParameterCount (type);
    22     if (!Nparams) {
    23         psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
    24         return(NULL);
    25     }     
    26 
    27     psf->params = psArrayAlloc (Nparams - 4);
    28     for (int i = 0; i < psf->params->n; i++) {
    29         // we need a way to set the allowed range of orders
    30         psf->params->data[i] = Polynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
    31     }
    32 
    33     psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
    34     return(psf);
    35 }
    36 
    37 static void pmPSF_TestFree (pmPSF_Test *test) {
     3// ********  pmPSFtry functions  **************************************************
     4// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
     5// * fitted version of the model, the pmPSF fit to the fitted model parameters,
     6// * and the PSF fits to the source. It also includes the statistics from the
     7// * fits, both the individual sources, and the collection
     8
     9// free a pmPSFtry structure
     10static void pmPSFtryFree (pmPSFtry *test) {
    3811
    3912  if (test == NULL) return;
    … …  
    4922}
    5023
    51 pmPSF_Test *pmPSF_TestAlloc (psArray *sources, char *modelName) {
    52 
    53     pmPSF_Test *test = (pmPSF_Test *) psAlloc(sizeof(pmPSF_Test));
     24// allocate a pmPSFtry based on the desired sources and the model (identified by name)
     25pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) {
     26
     27    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
    5428
    5529    // XXX probably need to increment ref counter
    … …  
    7448    }   
    7549
    76     psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree);
     50    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
    7751    return (test);
    7852}
    7953
    80 // test->mask values indicate reason source was rejected:
     54// build a pmPSFtry for the given model:
     55// - fit each source with the free-floating model
     56// - construct the pmPSF from the collection of models
     57// - fit each source with the PSF-parameter models
     58// - measure the pmPSF quality metric (dApResid)
     59
     60// sources used in for pmPSFtry may be masked by the analysis
     61// mask values indicate the reason the source was rejected:
    8162// 1: outlier in psf polynomial fit
    8263// 2: flt model failed to converge
    8364// 3: psf model failed to converge
    8465// 4: invalid source photometry
    85 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS)
     66// XXX EAM : use an enum for these values?
     67
     68pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS)
    8669{
    8770    bool status;
    … …  
    9376    int Npsf = 0;
    9477
    95     pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
     78    pmPSFtry *try = pmPSFtryAlloc (sources, modelName);
    9679
    9780    // stage 1:  fit an independent model (freeModel) to all sources
    9881    psTimerStart ("fit");
    99     for (int i = 0; i < test->sources->n; i++) {
    100 
    101         pmSource *source = test->sources->data[i];
    102         pmModel  *model  = pmSourceModelGuess (source, test->modelType);
     82    for (int i = 0; i < try->sources->n; i++) {
     83
     84        pmSource *source = try->sources->data[i];
     85        pmModel  *model  = pmSourceModelGuess (source, try->modelType);
    10386        x = source->peak->x;
    10487        y = source->peak->y;
    … …  
    11295        // exclude the poor fits
    11396        if (!status) {
    114           test->mask->data.U8[i] = 2;
     97          try->mask->data.U8[i] = 2;
    11598          psFree (model);
    11699          continue;
    117100        }
    118         test->modelFLT->data[i] = model;
     101        try->modelFLT->data[i] = model;
    119102        Nflt ++;
    120103    }
    121     psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    122     psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
     104    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     105    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
    123106
    124107    // make this optional?
    125     // DumpModelFits (test->modelFLT, "modelsFLT.dat");
     108    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
    126109
    127110    // stage 2: construct a psf (pmPSF) from this collection of model fits
    128     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
     111    pmPSFFromModels (try->psf, try->modelFLT, try->mask);
    129112
    130113    // stage 3: refit with fixed shape parameters
    131114    psTimerStart ("fit");
    132     for (int i = 0; i < test->sources->n; i++) {
     115    for (int i = 0; i < try->sources->n; i++) {
    133116        // masked for: bad model fit, outlier in parameters
    134         if (test->mask->data.U8[i]) continue;
    135 
    136         pmSource *source = test->sources->data[i];
    137         pmModel  *modelFLT = test->modelFLT->data[i];
     117        if (try->mask->data.U8[i]) continue;
     118
     119        pmSource *source = try->sources->data[i];
     120        pmModel  *modelFLT = try->modelFLT->data[i];
    138121
    139122        // set shape for this model based on PSF
    140         pmModel *modelPSF = pmModelFromPSF (modelFLT, test->psf);
     123        pmModel *modelPSF = pmModelFromPSF (modelFLT, try->psf);
    141124        x = source->peak->x;
    142125        y = source->peak->y;
    … …  
    147130        // skip poor fits
    148131        if (!status) {
    149             test->mask->data.U8[i] = 3;
     132            try->mask->data.U8[i] = 3;
    150133            psFree (modelPSF);
    151134            goto next_source;
    … …  
    153136
    154137        // otherwise, save the resulting model
    155         test->modelPSF->data[i] = modelPSF;
     138        try->modelPSF->data[i] = modelPSF;
    156139
    157140        // XXX : use a different aperture radius from the fit radius?
    … …  
    159142        // XXX : pass 'source' as input?
    160143        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
    161             test->mask->data.U8[i] = 4;
     144            try->mask->data.U8[i] = 4;
    162145            goto next_source;
    163146        }           
    164147
    165         test->metric->data.F64[i] = obsMag - fitMag;
    166         test->fitMag->data.F64[i] = fitMag;
     148        try->metric->data.F64[i] = obsMag - fitMag;
     149        try->fitMag->data.F64[i] = fitMag;
    167150        Npsf ++;
    168151
    … …  
    171154
    172155    }
    173     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
    174     psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
     156    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     157    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
    175158
    176159    // make this optional
    177     // DumpModelFits (test->modelPSF, "modelsPSF.dat");
    178 
    179     // XXX this function wants aperture radius from pmSourcePhotometry
    180     pmPSFMetricModel (test, RADIUS);
    181     psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n",
    182               modelName, test->ApResid, test->dApResid, test->skyBias);
    183 
    184     return (test);
    185 }
    186 
    187 // input: an array of pmModels, pre-allocated psf
    188 // some of the array entries may be NULL, ignore them
    189 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
    190 
    191     // construct the fit vectors from the collection of objects
    192     // use the mask to ignore missing fits
    193     psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
    194     psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
    195     psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
    196     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
    197 
    198     for (int i = 0; i < models->n; i++) {
    199         pmModel *model = models->data[i];
    200         if (model == NULL) continue;
    201 
    202         // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
    203         x->data.F64[i] = model->params->data.F32[2];
    204         y->data.F64[i] = model->params->data.F32[3];
    205     }
    206 
    207     // we are doing a robust fit.  after each pass, we drop points which are
    208     // more deviant than three sigma.
    209     // mask is currently updated for each pass. this is inconsistent?
    210 
    211     for (int i = 0; i < psf->params->n; i++) {
    212         for (int j = 0; j < models->n; j++) {
    213             pmModel *model = models->data[j];
    214             if (model == NULL) continue;
    215             z->data.F64[j] = model->params->data.F32[i + 4];
    216             dz->data.F64[j] = 1;
    217             // XXX EAM : need to use actual errors?
    218             // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
    219         }
    220 
    221         psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
    222         // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
    223         // psPolynomial2DDump (psf->params->data[i]);
    224     }
    225    
    226     psFree (x);
    227     psFree (y);
    228     psFree (z);
    229     psFree (dz);
    230     return (true);
    231 }
    232 
    233 pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) {
    234 
    235     // need to define the relationship between the modelFLT and modelPSF ?
    236    
    237     // find function used to set the model parameters
    238     pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
    239 
    240     // do we need a different model for the PSF vs FLT version?
    241     pmModel *modelPSF = pmModelAlloc (psf->type);
    242 
    243     // set model parameters for this source based on PSF information
    244     modelFromPSFFunc (modelPSF, modelFLT, psf);
    245 
    246     return (modelPSF);
    247 }
    248 
    249 bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
     160    // DumpModelFits (try->modelPSF, "modelsPSF.dat");
     161
     162    // XXX this function wants aperture radius for pmSourcePhotometry
     163    pmPSFtryMetric (try, RADIUS);
     164    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
     165              modelName, try->ApResid, try->dApResid, try->skyBias);
     166
     167    return (try);
     168}
     169
     170
     171bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
    250172
    251173  float dBin;
    252174  int   nKeep, nSkip;
    253175
    254   // the measured (aperture - fit) magnitudes (dA == test->metric)
     176  // the measured (aperture - fit) magnitudes (dA == try->metric)
    255177  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
    256178  //     dA = dAo + dsky/flux
    … …  
    262184
    263185  // rflux = ten(0.4*fitMag);
    264   psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
    265   for (int i = 0; i < test->sources->n; i++) {
    266     if (test->mask->data.U8[i]) continue;
    267     rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
     186  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     187  for (int i = 0; i < try->sources->n; i++) {
     188    if (try->mask->data.U8[i]) continue;
     189    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
    268190  }
    269191
    270192  // find min and max of (1/flux):
    271193  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    272   psVectorStats (stats, rflux, NULL, test->mask, 0xff);
     194  psVectorStats (stats, rflux, NULL, try->mask, 0xff);
    273195 
    274196  // build binned versions of rflux, metric
    … …  
    284206  for (int i = 0; i < daBin->n; i++) {
    285207
    286     psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
     208    psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
    287209    tmp->n = 0;
    288210
    289211    // accumulate data within bin range
    290     for (int j = 0; j < test->sources->n; j++) {
     212    for (int j = 0; j < try->sources->n; j++) {
    291213      // masked for: bad model fit, outlier in parameters
    292       if (test->mask->data.U8[j]) continue;
     214      if (try->mask->data.U8[j]) continue;
    293215   
    294216      // skip points with extreme dA values
    295       if (fabs(test->metric->data.F64[j]) > 0.5) continue;
     217      if (fabs(try->metric->data.F64[j]) > 0.5) continue;
    296218
    297219      // skip points outside of this bin
    … …  
    299221      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
    300222
    301       tmp->data.F64[tmp->n] = test->metric->data.F64[j];
     223      tmp->data.F64[tmp->n] = try->metric->data.F64[j];
    302224      tmp->n ++;
    303225    }
    … …  
    342264
    343265  // linear fit to rfBin, daBin
    344   psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
     266  psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
     267
     268  // XXX EAM : this is the intended API (cycle 7? cycle 8?)
     269  // poly = psVectorFitPolynomial1D (poly, maskB, 1, daBin, NULL, rfBin);
     270
     271  // XXX EAM : replace this when the above version is implemented
    345272  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
    346273
    347   psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
     274  psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
    348275  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
    349276
    … …  
    351278  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
    352279
    353   test->ApResid = poly->coeff[0];
    354   test->dApResid = stats->clippedStdev;
    355   test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
     280  try->ApResid = poly->coeff[0];
     281  try->dApResid = stats->clippedStdev;
     282  try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
    356283
    357284  psFree (rflux);
Note: See TracChangeset for help on using the changeset viewer.