IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26799


Ignore:
Timestamp:
Feb 5, 2010, 4:40:39 PM (16 years ago)
Author:
eugene
Message:

updates to follow changes to psphot APIs

Location:
branches/eam_branches/20091201/psastro/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psastro/src/psastroConvert.c

    r26259 r26799  
    11/** @file psastroConvert.c
    22 *
    3  *  @brief 
     3 *  @brief
    44 *
    55 *  @ingroup libpsastro
    … …  
    4747
    4848    // PSPHOT.SOURCES carries the pmSource objects (from psphot analysis or loaded externally)
    49     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    50     if (sources == NULL)
    51         return false;
     49    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     50    psAssert (detections, "missing detections?");
     51
     52    psArray *sources = detections->allSources;
     53    psAssert (sources, "missing sources?");
    5254
    5355    // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?)
    … …  
    6163    // XXX need to exit gracefully is inStars->n is 0 (or 1?)
    6264
    63     // we are going to select the brighter Nmax subset for astrometry 
     65    // we are going to select the brighter Nmax subset for astrometry
    6466    int nMax = psMetadataLookupS32 (&status, recipe, "PSASTRO.MAX.NRAW");
    6567    if (!status || !nMax) nMax = inStars->n;
    6668
    67     // we are going to select the brighter Nmax subset for astrometry 
     69    // we are going to select the brighter Nmax subset for astrometry
    6870    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.MAX.INST.MAG.RAW");
    6971    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.MIN.INST.MAG.RAW");
    7072
    71     // we are going to select the brighter Nmax subset for astrometry 
     73    // we are going to select the brighter Nmax subset for astrometry
    7274    pmSourceMode skip = PM_SOURCE_MODE_DEFAULT;
    7375    char *ignoreList = psMetadataLookupStr (&status, recipe, "PSASTRO.IGNORE");
    … …  
    7577      psArray *list = psStringSplitArray (ignoreList, ",", false);
    7678      for (int i = 0; i < list->n; i++) {
    77         pmSourceMode mode = pmSourceModeFromString (list->data[i]);
    78         if (mode == PM_SOURCE_MODE_DEFAULT) {
    79           psWarning ("unknown source mode in PSASTRO.IGNORE, skipping");
    80           continue;
    81         }
    82         skip |= mode;
     79        pmSourceMode mode = pmSourceModeFromString (list->data[i]);
     80        if (mode == PM_SOURCE_MODE_DEFAULT) {
     81          psWarning ("unknown source mode in PSASTRO.IGNORE, skipping");
     82          continue;
     83        }
     84        skip |= mode;
    8385      }
    8486      psFree (list);
    … …  
    9799
    98100    for (int i = 0; (i < inStars->n) && (j < rawStars->n); i++) {
    99         int n = index->data.S32[i];
    100         pmSource *source = sources->data[n];
     101        int n = index->data.S32[i];
     102        pmSource *source = sources->data[n];
    101103
    102104        psTrace ("psastro", 6, "mag: %f +/- %f, mode: %x, skip: %x\n", source->psfMag, source->errMag, source->mode, skip);
    103105
    104         if (source->mode & skip) {
    105           nModeSkip ++;
    106           continue;
    107         }
     106        if (source->mode & skip) {
     107          nModeSkip ++;
     108          continue;
     109        }
    108110
    109         if ((iMagMax != 0.0) && (source->psfMag > iMagMax)) {
    110           nFaintSkip ++;
    111           continue;
    112         }
    113         if ((iMagMin != 0.0) && (source->psfMag < iMagMin)) {
    114           nBrightSkip ++;
    115           continue;
    116         }
    117         if (!isfinite(source->psfMag)) {
    118           nInfSkip ++;
    119           continue;
    120         }
    121         mMin = PS_MIN (mMin, source->psfMag);
    122         mMax = PS_MAX (mMax, source->psfMag);
    123         rawStars->data[j] = psMemIncrRefCounter (inStars->data[n]);
    124         j++;
     111        if ((iMagMax != 0.0) && (source->psfMag > iMagMax)) {
     112          nFaintSkip ++;
     113          continue;
     114        }
     115        if ((iMagMin != 0.0) && (source->psfMag < iMagMin)) {
     116          nBrightSkip ++;
     117          continue;
     118        }
     119        if (!isfinite(source->psfMag)) {
     120          nInfSkip ++;
     121          continue;
     122        }
     123        mMin = PS_MIN (mMin, source->psfMag);
     124        mMax = PS_MAX (mMax, source->psfMag);
     125        rawStars->data[j] = psMemIncrRefCounter (inStars->data[n]);
     126        j++;
    125127    }
    126128    rawStars->n = j;
    … …  
    147149
    148150        // only accept the PSF sources?
    149         // XXX drop SATSTAR?
     151        // XXX drop SATSTAR?
    150152        // if (source->type != PM_SOURCE_TYPE_STAR) continue;
    151153
    … …  
    154156
    155157        psF32 *PAR = model->params->data.F32;
    156         psF32 *dPAR = model->dparams->data.F32;
     158        psF32 *dPAR = model->dparams->data.F32;
    157159
    158160        pmAstromObj *obj = pmAstromObjAlloc ();

  • branches/eam_branches/20091201/psastro/src/psastroMaskUpdates.c

    r26259 r26799  
    11/** @file psastroMaskUpdates.c
    22 *
    3  *  @brief 
     3 *  @brief
    44 *
    55 *  @ingroup libpsastro
    … …  
    3535    psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
    3636    psImageMaskType crosstalkMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for crosstalk ghosts
    37    
     37
    3838    // psImageMaskType maskBlank  = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels
    3939
    … …  
    5555      return(false);
    5656    }
    57    
     57
    5858    // convert star positions to ghost positions and add to the readout->analysis data
    5959    if (!psastroLoadGhosts (config)) {
    … …  
    8686    double REFSTAR_MASK_CROSSTALK_MAG_MAX  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX");
    8787    double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE");
    88    
     88
    8989    // select the input data sources
    9090    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    … …  
    188188                if (! readout->data_exists) { continue; }
    189189
    190                 // XXX why not do this?
    191                 pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
    192                 if (!readoutMask) continue;
     190                // XXX why not do this?
     191                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
     192                if (!readoutMask) continue;
    193193
    194194                // select the raw objects for this readout
    195                 // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
    196                 // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
    197                 // applied to the SUBSETs
     195                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
     196                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
     197                // applied to the SUBSETs
    198198                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    199199                if (!refstars) continue;
    … …  
    203203                // 2) diffraction spikes in direction ROT - ROTo
    204204                // 3) bleed trail in the direction of the readout
    205                 //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.
     205                //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.
    206206
    207207                for (int i = 0; i < refstars->n; i++) {
    208208                    pmAstromObj *ref = refstars->data[i];
    209                    
    210                     if (COUNT_GHOSTS) {
    211                         if (ref->Mag > GHOST_MAX_MAG) {
    212                             nGhosts ++;
    213                         }
    214                     }
     209
     210                    if (COUNT_GHOSTS) {
     211                        if (ref->Mag > GHOST_MAX_MAG) {
     212                            nGhosts ++;
     213                        }
     214                    }
    215215
    216216                    if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;
    217217
    218218
    219                     // the reference magnitudes have been converted from the instrumental
    220                     // values supplied in the recipe to apparent mags (above)
     219                    // the reference magnitudes have been converted from the instrumental
     220                    // values supplied in the recipe to apparent mags (above)
    221221
    222222                    // CIRCLE around the stars (scaled by magnitude)
    … …  
    224224
    225225                    // XXX for now, assume cell binning is 1x1 relative to chip
    226                     psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
     226                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
    227227
    228228                    for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
    … …  
    230230                        float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
    231231
    232                         // LINE for boundaries of the saturation spikes (scaled by magnitude)
    233                         // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
    234 
    235                         float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
    236                         float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
    237                         float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
    238                         // XXX we can make the width depend on the spike as well...
    239                         // The length should also be a function of the image background level
    240 
    241                         psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
    242                     }
    243 
    244                     // This masking option was needed for persistent charge trails in GPC1; it
    245                     // has since been replaced with 'burntool', which is applied upon readout
    246                     // by the camera software, and therefore is aware of the image sequence.
    247                     if (REFSTAR_MASK_BLEED) {
    248                         // convert x,y chip coordinates to cells in maskChip
    249                         pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
    250 
    251                         // LINE for boundaries of the bleed lines
    252                         if (refCell) {
    253                             float xCell = 0.0;
    254                             float yCell = 0.0;
    255                             float xEnd = 0.0;
    256                             float yEnd = 0.0;
    257                             // find coordinate of star on cell
    258                             pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
    259                             // find coordinate of end-point on chip
    260 
    261                             int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
    262                             pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
    263 
    264                             float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
    265                             psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
    266                         }
    267                     }
     232                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
     233                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
     234
     235                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
     236                        float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
     237                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
     238                        // XXX we can make the width depend on the spike as well...
     239                        // The length should also be a function of the image background level
     240
     241                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
     242                    }
     243
     244                    // This masking option was needed for persistent charge trails in GPC1; it
     245                    // has since been replaced with 'burntool', which is applied upon readout
     246                    // by the camera software, and therefore is aware of the image sequence.
     247                    if (REFSTAR_MASK_BLEED) {
     248                        // convert x,y chip coordinates to cells in maskChip
     249                        pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
     250
     251                        // LINE for boundaries of the bleed lines
     252                        if (refCell) {
     253                            float xCell = 0.0;
     254                            float yCell = 0.0;
     255                            float xEnd = 0.0;
     256                            float yEnd = 0.0;
     257                            // find coordinate of star on cell
     258                            pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
     259                            // find coordinate of end-point on chip
     260
     261                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
     262                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
     263
     264                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
     265                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
     266                        }
     267                    }
    268268                }
    269269
    … …  
    274274                // supplied here.
    275275                psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
    276                 if (ghosts) { 
    277                     // mask the ghosts on this readout
    278                     for (int i = 0; i < ghosts->n; i++) {
    279                         psastroGhost *ghost = ghosts->data[i];
    280                         // XXX bright vs faint ghost bits? (OR with SUSPECT)
    281                         psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
    282                     }
    283                 }
     276                if (ghosts) {
     277                    // mask the ghosts on this readout
     278                    for (int i = 0; i < ghosts->n; i++) {
     279                        psastroGhost *ghost = ghosts->data[i];
     280                        // XXX bright vs faint ghost bits? (OR with SUSPECT)
     281                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
     282                    }
     283                }
    284284
    285285                // Select the glint mask regions for this readout (loaded in
    … …  
    290290                psArray *glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS");
    291291                if (glints) {
    292                     // mask the glints on this readout
    293                     for (int i = 0; i < glints->n; i++) {
    294                         psRegion *glint = glints->data[i];
    295                         psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
    296                     }
    297                 }
    298 
    299                 // Select the crosstalk artifact regions for this readout, and mask a circular region
    300                 // corresponding to the source star
    301 
    302                 psArray *crosstalks = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.CROSSTALKS");
    303                 if (crosstalks) {
    304                   for (int i = 0; i < crosstalks->n; i++) {
    305                     pmAstromObj *ref = crosstalks->data[i];
    306                     float radius = REFSTAR_MASK_CROSSTALK_MAG_SLOPE * (REFSTAR_MASK_CROSSTALK_MAG_MAX - ref->Mag);
    307                     psTrace("psastro.crosstalk",2,"Masking star on Chip %s @ (%f,%f) Magnitude: %f Radius %f\n",
    308                             psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"),
    309                             ref->chip->x,ref->chip->y,ref->Mag,radius);
    310                     // XXX for now, assume cell binning is 1x1 relative to chip
    311                     psastroMaskCircle (readoutMask->mask, crosstalkMaskValue, ref->chip->x, ref->chip->y, radius, radius);
    312                   }
    313                 }
    314                
    315                 // this probably should move into a function of its own:
    316                 {
    317                     // select the raw objects for this readout, flag is they fall in a mask
    318                     psArray *inSources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    319                     if (inSources == NULL) continue;
    320                
    321                     // create a replacement output array:
    322                     // psArray *outSources = psAllocArrayEmpty(100);
    323 
    324                     // XXX finish this: raise a bit for stars that land on certain types of masks;
    325                     // others (eg, bright star core) should be ignored.
    326                     for (int i = 0; i < inSources->n; i++) {
    327                         pmSource *source = inSources->data[i];
    328 
    329                         int xChip = source->peak->x;
    330                         int yChip = source->peak->y;
    331 
    332                         bool onChip = true;
    333                         onChip &= (xChip >= 0);
    334                         onChip &= (xChip < readoutMask->mask->numCols);
    335                         onChip &= (yChip >= 0);
    336                         onChip &= (yChip < readoutMask->mask->numRows);
    337                         if (!onChip) {
    338                             // if the source is off the edge of the chip, raise a different bit?
    339                             source->mode |= PM_SOURCE_MODE_OFF_CHIP;
    340                             continue;
    341                         }
    342 
    343                         psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
    344                         if (value & ghostMaskValue) {
    345                             source->mode |= PM_SOURCE_MODE_ON_GHOST;
    346                         }
    347                         // XXX note that for now, glint and ghost are identical
    348                         pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
    349                         if (value & glintMaskValue) {
    350                             source->mode |= PM_SOURCE_MODE_ON_GLINT;
    351                         }
    352                         if (value & spikeMaskValue) {
    353                             source->mode |= PM_SOURCE_MODE_ON_SPIKE;
    354                         }
    355                     }
    356                 }
     292                    // mask the glints on this readout
     293                    for (int i = 0; i < glints->n; i++) {
     294                        psRegion *glint = glints->data[i];
     295                        psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
     296                    }
     297                }
     298
     299                // Select the crosstalk artifact regions for this readout, and mask a circular region
     300                // corresponding to the source star
     301
     302                psArray *crosstalks = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.CROSSTALKS");
     303                if (crosstalks) {
     304                  for (int i = 0; i < crosstalks->n; i++) {
     305                    pmAstromObj *ref = crosstalks->data[i];
     306                    float radius = REFSTAR_MASK_CROSSTALK_MAG_SLOPE * (REFSTAR_MASK_CROSSTALK_MAG_MAX - ref->Mag);
     307                    psTrace("psastro.crosstalk",2,"Masking star on Chip %s @ (%f,%f) Magnitude: %f Radius %f\n",
     308                            psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"),
     309                            ref->chip->x,ref->chip->y,ref->Mag,radius);
     310                    // XXX for now, assume cell binning is 1x1 relative to chip
     311                    psastroMaskCircle (readoutMask->mask, crosstalkMaskValue, ref->chip->x, ref->chip->y, radius, radius);
     312                  }
     313                }
     314
     315                // this probably should move into a function of its own:
     316                {
     317                    // select the raw objects for this readout, flag is they fall in a mask
     318                    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     319                    if (!detections) continue;
     320
     321                    psArray *inSources = detections->allSources;
     322                    psAssert (inSources, "missing sources?");
     323
     324                    // create a replacement output array:
     325                    // psArray *outSources = psAllocArrayEmpty(100);
     326
     327                    // XXX finish this: raise a bit for stars that land on certain types of masks;
     328                    // others (eg, bright star core) should be ignored.
     329                    for (int i = 0; i < inSources->n; i++) {
     330                        pmSource *source = inSources->data[i];
     331
     332                        int xChip = source->peak->x;
     333                        int yChip = source->peak->y;
     334
     335                        bool onChip = true;
     336                        onChip &= (xChip >= 0);
     337                        onChip &= (xChip < readoutMask->mask->numCols);
     338                        onChip &= (yChip >= 0);
     339                        onChip &= (yChip < readoutMask->mask->numRows);
     340                        if (!onChip) {
     341                            // if the source is off the edge of the chip, raise a different bit?
     342                            source->mode |= PM_SOURCE_MODE_OFF_CHIP;
     343                            continue;
     344                        }
     345
     346                        psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
     347                        if (value & ghostMaskValue) {
     348                            source->mode |= PM_SOURCE_MODE_ON_GHOST;
     349                        }
     350                        // XXX note that for now, glint and ghost are identical
     351                        pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
     352                        if (value & glintMaskValue) {
     353                            source->mode |= PM_SOURCE_MODE_ON_GLINT;
     354                        }
     355                        if (value & spikeMaskValue) {
     356                            source->mode |= PM_SOURCE_MODE_ON_SPIKE;
     357                        }
     358                    }
     359                }
    357360            }
    358361        }
    … …  
    375378
    376379    if (COUNT_GHOSTS) {
    377         // save nGhosts to update header.
    378         psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
    379         if (!updates) {
    380             updates = psMetadataAlloc ();
    381             psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
    382             psFree (updates);
    383         }
    384         psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
     380        // save nGhosts to update header.
     381        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
     382        if (!updates) {
     383            updates = psMetadataAlloc ();
     384            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
     385            psFree (updates);
     386        }
     387        psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
    385388    }
    386389
Note: See TracChangeset for help on using the changeset viewer.