IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5654


Ignore:
Timestamp:
Dec 1, 2005, 9:55:16 AM (20 years ago)
Author:
eugene
Message:

adding ensemble photometry

Location:
trunk/psphot/src
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psSparse.c

    r5653 r5654  
    11# include "psphot.h"
    2 
    3 typedef struct {
    4     psVector *Aij;
    5     psVector *Bfj;
    6     psVector *Qii;
    7     psVector *Si;
    8     psVector *Sj;
    9     int Nelem;
    10     int Nrows;
    11 } psSparse;
    122
    133void psSparseMatrixTest () {
    … …  
    166    psSparse *sparse = psSparseAlloc (3, 9);
    177
    18     psSparseMatrixElement (sparse, 0, 0, 1.0);
    19     psSparseMatrixElement (sparse, 1, 1, 1.0);
     8    psSparseMatrixElement (sparse, 0, 0, 3.0);
     9    psSparseMatrixElement (sparse, 1, 1, 2.0);
    2010    psSparseMatrixElement (sparse, 2, 2, 1.0);
    2111
    22     psSparseMatrixElement (sparse, 0, 1, 0.1);
    23     psSparseMatrixElement (sparse, 0, 2, 0.1);
    24 
    25     psSparseMatrixResort (sparse);
     12    psSparseMatrixElement (sparse, 1, 0, 0.1);
     13    psSparseMatrixElement (sparse, 2, 0, -0.1);
     14
     15    psSparseResort (sparse);
    2616    for (int i = 0; i < sparse->Nelem; i++) {
    2717        fprintf (stderr, "%d %d %f\n",
    … …  
    3323    psVector *x = psVectorAlloc (3, PS_DATA_F32);
    3424    x->data.F32[0] = 3;
    35     x->data.F32[0] = 5;
    36     x->data.F32[0] = 7;
    37 
    38     psVector B = psSparseMatrixTimesVector (NULL, sparse, x);
     25    x->data.F32[1] = 5;
     26    x->data.F32[2] = 7;
     27    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     28
     29    psVector *B = psSparseMatrixTimesVector (NULL, sparse, x);
    3930    fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]);
    4031
    … …  
    4334    sparse->Bfj->data.F32[2] = B->data.F32[2];
    4435
    45     x = psSparseSolve (x, sparse);
     36    x = psSparseSolve (x, sparse, 0);
     37    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     38
     39    x = psSparseSolve (x, sparse, 1);
     40    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     41
     42    x = psSparseSolve (x, sparse, 2);
     43    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     44
     45    x = psSparseSolve (x, sparse, 3);
     46    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     47
     48    x = psSparseSolve (x, sparse, 4);
    4649    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    4750    return;
    … …  
    5053void psSparseResort (psSparse *sparse) {
    5154
    52     Nelem = sparse->Nelem;
     55    int Nelem = sparse->Nelem;
    5356
    5457    psVector *index = psVectorSortIndex (NULL, sparse->Sj);
    55     psVector Aij = sparse->Aij;
    56     psVector Si = sparse->Si;
    57     psVector Sj = sparse->Sj;
     58    psVector *Aij = sparse->Aij;
     59    psVector *Si = sparse->Si;
     60    psVector *Sj = sparse->Sj;
    5861
    5962    // allocate new temporary vectors
    … …  
    6164    psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
    6265    psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
    63     for (i = 0; i < Nelem; i++) {
    64         j = index->data.U32[i];
     66    for (int i = 0; i < Nelem; i++) {
     67        int j = index->data.U32[i];
    6568        tAij->data.F32[i] = Aij->data.F32[j];
    6669        tSi->data.S32[i]  = Si->data.S32[j];
    … …  
    8083void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) {
    8184
     85    int k;
     86
    8287    if (i < j) {
    8388        fprintf (stderr, "*** error: subdiagonal element ***\n");
    … …  
    9196        // check vectors lengths and extend if needed
    9297        if (sparse->Nelem >= sparse->Aij->nalloc) {
    93             sparse->Aij->nalloc += 100;     
    94             psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
    95             sparse->Si->nalloc += 100;     
    96             psVectorRealloc (sparse->Si, sparse->Si->nalloc);
    97             sparse->Sj->nalloc += 100;     
    98             psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
     98            psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
     99            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
     100            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
    99101        }
    100102
    … …  
    111113        // check vectors lengths and extend if needed
    112114        if (sparse->Nelem >= sparse->Aij->nalloc - 1) {
    113             sparse->Aij->nalloc += 100;     
    114             psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
    115             sparse->Si->nalloc += 100;     
    116             psVectorRealloc (sparse->Si, sparse->Si->nalloc);
    117             sparse->Sj->nalloc += 100;     
    118             psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
     115            psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
     116            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
     117            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
    119118        }
    120119
    … …  
    139138void psSparseVectorElement (psSparse *sparse, int i, float value) {
    140139
    141     sparse->Bjf->data.F32[i] = value;
     140    sparse->Bfj->data.F32[i] = value;
    142141    return;
    143142}
    … …  
    146145psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) {
    147146
     147    int i, Nelem;
     148    float F;
     149
    148150    if (output == NULL) {
    149         psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
     151        output = psVectorAlloc (vector->n, PS_DATA_F32);
    150152    }
    151153
    … …  
    154156        F = 0;
    155157        while (matrix->Sj->data.S32[Nelem] == j) {
    156             i = matrix->Sj->data.S32[Nelem];
    157             F += vector->data.F32[i] * matrix->Aij->data.F32[i];
     158            i = matrix->Si->data.S32[Nelem];
     159            F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem];
     160            Nelem++;
    158161        }
    159162        output->data.F32[j] = F;
    … …  
    162165}
    163166
    164 psVector *psSparseSolve (psVector *guess, psSparse *sparse) {
     167psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
    165168
    166169    psF32 dG;
    … …  
    169172    psVector *Bfj = sparse->Bfj;
    170173
    171     guess = psVectorCopy (guess, Bfj);
     174    guess = psVectorCopy (guess, Bfj, PS_DATA_F32);
    172175
    173176    // temporary storage for intermediate results
    174177    psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
    175178
    176     for (int j = 0; j < 2; j++) {
     179    for (int j = 0; j < Niter; j++) {
    177180        dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
    178         for (int i = 0; i < dG->n; i++) {
     181        for (int i = 0; i < dQ->n; i++) {
    179182            dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
    180183            guess->data.F32[i] -= dG;
    … …  
    187190static void psSparseFree (psSparse *sparse) {
    188191    if (sparse == NULL) return;
    189     psFree (sparse.Aij);
    190     psFree (sparse.Bfj);
    191     psFree (sparse.Qii);
    192     psFree (sparse.Si);
    193     psFree (sparse.Sj);
     192    psFree (sparse->Aij);
     193    psFree (sparse->Bfj);
     194    psFree (sparse->Qii);
     195    psFree (sparse->Si);
     196    psFree (sparse->Sj);
    194197    return;
    195198}
    … …  
    212215
    213216    sparse->Nrows = Nrows;
     217    psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree);
    214218    return (sparse);
    215219}

  • trunk/psphot/src/psphot.h

    r5643 r5654  
    99# include "pmPSFtry.h"
    1010# include "pmModelGroup.h"
     11# include "psSparse.h"
    1112
    1213typedef struct {

  • trunk/psphot/src/psphotSparseMatrix.c

    r5653 r5654  
    11# include "psphot.h"
    22
    3 psSparse *psphotStarOverlaps (psArray *sources) {
     3bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) {
    44
    5     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    6     psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    7     psVector *params = model->params;
     5    bool  status;
     6    float x;
     7    float y;
     8    float Sky;
     9    int   Nfit = 0;
    810
    9     // source analysis is done in S/N order (brightest first)
     11    psTimerStart ("psphot");
     12
     13    // source analysis is done in spatial order
    1014    sources = psArraySort (sources, psphotSortByY);
    1115
    12     // fill out the models for each object
    13     // XXX need to assign the model parameters using the PSF function, with unit peak
     16    float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     17    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
     18    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
     19
     20    // set the object surface-brightness limit for fitted pixels
     21    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
     22    psLogMsg ("psphot.apply_psf_model", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
     23
     24    // vector for model input coordinates
     25    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     26
     27    // this function specifies the radius at this the model hits the given flux
     28    pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type);
     29
     30    // this is the function to get the model function
     31    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     32
     33    // pre-calculate all model pixels
    1434    psArray  *pixlist = psArrayAlloc (sources->n);
    1535    for (int k = 0; k < sources->n; k++) {
    16         src = sources->data[k];
     36        source = sources->data[k];
     37
     38        // use the source moments, etc to guess basic model parameters
     39        pmModel *modelFLT = pmSourceModelGuess (source, psf->type);
     40
     41        // set PSF parameters for this model
     42        pmModel *model = pmModelFromPSF (modelFLT, psf);
     43        // model->params->data.F32[2] = source->peak->x;
     44        // model->params->data.F32[3] = source->peak->y;
     45        model->params->data.F32[1] = 1.0;
     46        x = model->params->data.F32[2];
     47        y = model->params->data.F32[3];
     48        psFree (modelFLT);
     49
     50        // set the fit radius based on the object flux limit and the model
     51        // XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
     52        model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
     53        if (isnan(model->radius)) {
     54          psAbort ("apply_psf_model", "error in radius");
     55        }
     56
     57        // build the model image
    1758        nx = src->pixels->numCols;
    1859        ny = src->pixels->numRows;
    1960        im = psImageAlloc (nx, ny, PS_DATA_F32);
     61        mk = psImageAlloc (nx, ny, PS_DATA_U8);
     62
     63        // fill in the model pixel values
     64        psImageKeepCircle (mk, x, y, model->radius, "OR", PSPHOT_MASK_MARKED);
     65        psVector *params = model->params;
    2066        for (int i = 0; i < nx; i++) {
    2167            for (int j = 0; j < ny; j++) {
     68                if (mk->data.U8[j][i]) continue;
    2269                x->data.F32[0] = (float) i + pix->col0;
    2370                x->data.F32[1] = (float) j + pix->row0;
    … …  
    2875    }
    2976   
    30     psVector *Aij = psVectorAlloc (100, PS_DATA_F32);
    31     psVector *Bfj = psVectorAlloc (100, PS_DATA_F32);
    32     psVector *Si  = psVectorAlloc (100, PS_DATA_S32);
    33     psVector *Sj  = psVectorAlloc (100, PS_DATA_S32);
     77    // fill out the sparse matrix
     78    psSparse *sparse = psSparseAlloc (sources->n, 100);
     79    psVector *Aij = sparse->Aij;
     80    psVector *Bfj = sparse->Bfj;
     81    psVector *Si  = sparse->Si;
     82    psVector *Sj  = sparse->Sj;
    3483
    3584    int Ri, Rj;
    … …  
    4291        // diagonal element (auto-cross-product)
    4392        f = psphotCrossProduct (pi, pi);
    44         Aij->data.F32[Nelem] = f;
    45         Si->data.S32[Nelem] = i;
    46         Sj->data.S32[Nelem] = j;
    47         Nelem ++;
    48 
    49         Qii->data.F32[i] = f;
     93        psSparseMatrixElement (sparse, i, i, f);
    5094
    5195        // find the image x model value
    5296        f = psphotCrossProduct (pi, pf);
    53         Bfj->data.F32[i] = f;
     97        psSparseVectorElement (sparse, i, f);
    5498
    5599        // loop over all other stars following this one
    … …  
    64108            // got an overlap; calculate cross-product and add to output array
    65109            f = psphotCrossProduct (pi, pj);
    66            
    67             Aij->data.F32[Nelem] = f;
    68             Si->data.S32[Nelem] = i;
    69             Sj->data.S32[Nelem] = j;
    70             Nelem ++;
    71 
    72             Aij->data.F32[Nelem] = f;
    73             Si->data.S32[Nelem] = j;
    74             Sj->data.S32[Nelem] = i;
    75             Nelem ++;
     110            psSparseMatrixElement (sparse, j, i, f);
    76111        }
    77112    }
    78 
    79     index = psVectorSortIndex (NULL, Sj);
    80 
    81     psVector *tAij = psVectorAlloc (Nelem, PS_DATA_F32);
    82     psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
    83     psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
    84     for (i = 0; i < Nelem; i++) {
    85         j = index->data.U32[i];
    86         tAij->data.F32[i] = Aij->data.F32[j];
    87         tSi->data.S32[i]  = Si->data.S32[j];
    88         tSj->data.S32[i]  = Sj->data.S32[j];
    89     }
    90     psFree (Aij);
    91     psFree (Si);
    92     psFree (Sj);
    93 
    94     sparse->Aij = Aij;
    95     sparse->Bfj = Bfj;
    96     sparse->Qii = Qii;
    97     sparse->Si = Si;
    98     sparse->Sj = Sj;
    99     sparse->Nelem = Nelem;
    100     sparse->Nrows = sources->n;
    101 
     113   
     114    psSparseResort (sparse);
    102115    return (match);
    103116}

  • trunk/psphot/src/psphotTest.c

    r5652 r5654  
    33int main (int argc, char **argv) {
    44
    5   psphotSparseMatrixTest ();
     5  psSparseMatrixTest ();
    66
    77}
Note: See TracChangeset for help on using the changeset viewer.