IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 11, 2006, 5:27:13 PM (20 years ago)
Author:
eugene
Message:

major revisions to use the new pmFPAfile,view I/O handling

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psModulesUtils.c

    r6522 r6571  
    11# include "psModulesUtils.h"
    2 
    3 // extract config informatin from config data or from image header, if specified
    4 // XXX EAM : this function should be replaced with Paul's image/header/metadata load scheme
    5 psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
    6 
    7     char *source;
    8     char *keyword;
    9     psF32 value;
    10     psMetadataItem *item;
    11 
    12     // look for the entry in the config collection
    13     item = psMetadataLookup (config, name);
    14     if (item == NULL) {
    15         psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
    16         value = psMetadataLookupF32 (status, header, name);
    17         if (*status == false) {
    18             psAbort ("pmConfigLookupF32", "no key %s in header", name);
    19         }
    20         *status = true;
    21         return (value);
    22     }   
    23 
    24     // I'm either expecting a string, with the name "HD:keyword"...
    25     if (item->type == PS_DATA_STRING) {
    26         source = item->data.V;
    27         if (!strncasecmp (source, "HD:", 3)) {
    28             keyword = &source[3];
    29             value = psMetadataLookupF32 (status, header, keyword);
    30             if (*status == false) {
    31                 psAbort ("pmConfigLookupF32", "no key %s in config", name);
    32             }
    33             *status = true;
    34             // psFree (item);
    35             return (value);
    36         }       
    37     }
    38 
    39     //  ... or a value (F32?)
    40     if (item->type == PS_DATA_F32) {
    41         value = item->data.F32;
    42         // psFree (item);
    43         return (value);
    44     }
    45 
    46     *status = false;
    47     psError(PS_ERR_UNKNOWN, true, "invalid item");
    48     return (0);
    49 }
    50 
    51 bool pmModelFitStatus (pmModel *model) {
    52 
    53     bool status;
    54 
    55     pmModelFitStatusFunc statusFunc = pmModelFitStatusFunc_GetFunction (model->type);
    56     status = statusFunc (model);
    57 
    58     return (status);
    59 }
    60 
    61 pmModel *pmModelSelect (pmSource *source) {
    62     switch (source->type) {
    63       case PM_SOURCE_STAR:
    64         return source->modelPSF;
    65        
    66       case PM_SOURCE_EXTENDED:
    67         return source->modelEXT;
    68        
    69       default:
    70         return NULL;
    71     }
    72     return NULL;
    73 }
    74 
    75 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {
    76 
    77     float obsSum = 0;
    78     float fitSum = 0;
    79     float sky = model->params->data.F32[0];
    80 
    81     *fitMag = 99.0;
    82     *obsMag = 99.0;
    83 
    84     pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    85     fitSum = modelFluxFunc (model->params);
    86     if (fitSum <= 0) return false;
    87     if (!isfinite(fitSum)) return false;
    88     *fitMag = -2.5*log10(fitSum);
    89 
    90     for (int ix = 0; ix < image->numCols; ix++) {
    91         for (int iy = 0; iy < image->numRows; iy++) {
    92             if (mask->data.U8[iy][ix]) continue;
    93             obsSum += image->data.F32[iy][ix] - sky;
    94         }
    95     }
    96     if (obsSum <= 0) return false;
    97     *obsMag = -2.5*log10(obsSum);
    98 
    99     return (true);
    100 }
    101 
    102 // force the fitted sky to meet the source flux at outer radius
    103 bool pmModelSkyOffset (pmModel *model, float x, float y, float radius) {
    104 
    105     float flux;
    106 
    107     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    108 
    109     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    110     psVector *params = model->params;
    111  
    112     coord->data.F32[0] = x + radius;
    113     coord->data.F32[1] = y;
    114     flux = modelFunc (NULL, params, coord);
    115     params->data.F32[0] = flux;
    116 
    117     psFree (coord);
    118 
    119     return true;
    120 }
    121 
    122 pmModel *pmModelCopy (pmModel *model) {
    123 
    124     pmModel *new = pmModelAlloc (model->type);
    125    
    126     new->chisq  = model->chisq;
    127     new->nDOF   = model->nDOF;
    128     new->nIter  = model->nIter;
    129     new->status = model->status;
    130     new->radius = model->radius;
    131 
    132     for (int i = 0; i < new->params->n; i++) {
    133         new->params->data.F32[i]  = model->params->data.F32[i];
    134         new->dparams->data.F32[i] = model->dparams->data.F32[i];
    135     }
    136    
    137     return (new);
    138 }
    139 
    140 float pmSourceCrossProduct (pmSource *Mi, pmSource *Mj) {
    141 
    142     int Xs, Xe, Ys, Ye;
    143     int xi, xj, yi, yj;
    144     int xIs, xJs, yIs, yJs;
    145     int xIe, yIe;
    146     float flux, wt;
    147 
    148     psImage *Pi = Mi->pixels;
    149     psImage *Pj = Mj->pixels;
    150 
    151     psImage *Wi = Mi->weight;
    152 
    153     psImage *Ti = Mi->mask;
    154     psImage *Tj = Mj->mask;
    155 
    156     Xs = PS_MAX (Pi->col0, Pj->col0);
    157     Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
    158    
    159     Ys = PS_MAX (Pi->row0, Pj->row0);
    160     Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
    161    
    162     xIs = Xs - Pi->col0;
    163     xJs = Xs - Pj->col0;
    164     yIs = Ys - Pi->row0;
    165     yJs = Ys - Pj->row0;
    166 
    167     xIe = Xe - Pi->col0;
    168     yIe = Ye - Pi->row0;
    169 
    170     // note that this is addressing the same image pixels,
    171     // though only if both are source not model images
    172     flux = 0;
    173     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    174         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    175             if (Ti->data.U8[yi][xi]) continue;
    176             if (Tj->data.U8[yj][xj]) continue;
    177             wt = Wi->data.F32[yi][xi];
    178             if (wt > 0) {
    179                 flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
    180             }
    181         }
    182     }
    183     return (flux);
    184 }
    185 
    186 float pmSourceCrossWeight (pmSource *Mi, pmSource *Mj) {
    187 
    188     int Xs, Xe, Ys, Ye;
    189     int xi, xj, yi, yj;
    190     int xIs, xJs, yIs, yJs;
    191     int xIe, yIe;
    192     float flux, wt;
    193 
    194     psImage *Pi = Mi->pixels;
    195     psImage *Pj = Mj->pixels;
    196 
    197     psImage *Wi = Mi->weight;
    198 
    199     psImage *Ti = Mi->mask;
    200     psImage *Tj = Mj->mask;
    201 
    202     Xs = PS_MAX (Pi->col0, Pj->col0);
    203     Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
    204    
    205     Ys = PS_MAX (Pi->row0, Pj->row0);
    206     Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
    207    
    208     xIs = Xs - Pi->col0;
    209     xJs = Xs - Pj->col0;
    210     yIs = Ys - Pi->row0;
    211     yJs = Ys - Pj->row0;
    212 
    213     xIe = Xe - Pi->col0;
    214     yIe = Ye - Pi->row0;
    215 
    216     // note that this is addressing the same image pixels,
    217     // though only if both are source not model images
    218     flux = 0;
    219     for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
    220         for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
    221             if (Ti->data.U8[yi][xi]) continue;
    222             if (Tj->data.U8[yj][xj]) continue;
    223             wt = Wi->data.F32[yi][xi];
    224             if (wt > 0) {
    225                 flux += 1.0 / wt;
    226             }
    227         }
    228     }
    229     return (flux);
    230 }
    2312
    2323static void ppConfigFree (ppConfig *config) {
Note: See TracChangeset for help on using the changeset viewer.