Changeset 27838 for branches/tap_branches/psastro/src/psastroMaskUpdates.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psastro/src/psastroMaskUpdates.c (modified) (11 diffs)

branches/tap_branches/psastro/src/psastroMaskUpdates.c

@@ -1 +1 @@
 /** @file psastroMaskUpdates.c
  *
- *  @brief 
+ *  @brief
  *
  *  @ingroup libpsastro
@@ -34 +34 @@
     psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
     psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
+    psImageMaskType crosstalkMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for crosstalk ghosts
 
     // psImageMaskType maskBlank  = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels
@@ -47 +48 @@
     bool REFSTAR_MASK                      = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK");
     if (!REFSTAR_MASK) return true;
+
+    // convert star positions to crosstalk artifact positions and add to the readout->analysis data
+    bool CROSSTALK_MASK                    = psMetadataLookupBool (&status, recipe, "CROSSTALK_MASK");
+    if (CROSSTALK_MASK && !psastroLoadCrosstalk (config)) {
+      psError(PSASTRO_ERR_CONFIG, false, "Error loading crosstalk data");
+      return(false);
+    }
 
     // convert star positions to ghost positions and add to the readout->analysis data
@@ -73 +81 @@
     double REFSTAR_MASK_SATSPIKE_MAG_MAX   = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_MAG_MAX");
     double REFSTAR_MASK_SATSPIKE_WIDTH     = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_WIDTH");
+    double REFSTAR_MASK_SATSPIKE_OFFSET    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSPIKE_OFFSET");
     double REFSTAR_MASK_BLEED_MAG_MAX      = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX");
     double REFSTAR_MASK_BLEED_MAG_SLOPE    = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE");
+
+    double REFSTAR_MASK_CROSSTALK_MAG_MAX  = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX");
+    double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE");
 
     // select the input data sources
@@ -177 +189 @@
                 if (! readout->data_exists) { continue; }
 
-                // XXX why not do this? 
-                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
-                if (!readoutMask) continue;
+                // XXX why not do this?
+                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
+                if (!readoutMask) continue;
 
                 // select the raw objects for this readout
+                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
+                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
+                // applied to the SUBSETs
                 psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
                 if (!refstars) continue;
@@ -189 +204 @@
                 // 2) diffraction spikes in direction ROT - ROTo
                 // 3) bleed trail in the direction of the readout
-                //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.
+                //    update: with 'burntool' applied to the data, the bleed trail mask is not needed.
 
                 for (int i = 0; i < refstars->n; i++) {
                     pmAstromObj *ref = refstars->data[i];
-                    
-                    if (COUNT_GHOSTS) {
-                        if (ref->Mag > GHOST_MAX_MAG) {
-                            nGhosts ++;
-                        }
-                    }
+
+                    if (COUNT_GHOSTS) {
+                        if (ref->Mag > GHOST_MAX_MAG) {
+                            nGhosts ++;
+                        }
+                    }
 
                     if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue;
 
-                    // the reference magnitudes have been converted from the instrumental
-                    // values supplied in the recipe to apparent mags (above)
+
+                    // the reference magnitudes have been converted from the instrumental
+                    // values supplied in the recipe to apparent mags (above)
 
                     // CIRCLE around the stars (scaled by magnitude)
@@ -209 +225 @@
 
                     // XXX for now, assume cell binning is 1x1 relative to chip
-                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
+                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
 
                     for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
@@ -215 +231 @@
                         float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
 
-                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
-                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
-
-                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
-                        float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
-                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
-                        // XXX we can make the width depend on the spike as well...
-                        // The length should also be a function of the image background level
-
-                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
-                    }
-
-                    // This masking option was needed for persistent charge trails in GPC1; it
-                    // has since been replaced with 'burntool', which is applied upon readout
-                    // by the camera software, and therefore is aware of the image sequence.
-                    if (REFSTAR_MASK_BLEED) {
-                        // convert x,y chip coordinates to cells in maskChip
-                        pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
-
-                        // LINE for boundaries of the bleed lines
-                        if (refCell) {
-                            float xCell = 0.0;
-                            float yCell = 0.0;
-                            float xEnd = 0.0;
-                            float yEnd = 0.0;
-                            // find coordinate of star on cell
-                            pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
-                            // find coordinate of end-point on chip
-
-                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
-                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
-
-                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
-                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
+                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
+                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
+
+                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
+                        float spikeLength = pow(10,REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag)) - REFSTAR_MASK_SATSPIKE_OFFSET;
+                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
+                        if (spikeLength < 0.0) {
+                          spikeLength = 0.0;
                         }
-                    }
+                        // XXX we can make the width depend on the spike as well...
+                        // The length should also be a function of the image background level
+
+                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
+                    }
+
+                    // This masking option was needed for persistent charge trails in GPC1; it
+                    // has since been replaced with 'burntool', which is applied upon readout
+                    // by the camera software, and therefore is aware of the image sequence.
+                    if (REFSTAR_MASK_BLEED) {
+                        // convert x,y chip coordinates to cells in maskChip
+                        pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y);
+
+                        // LINE for boundaries of the bleed lines
+                        if (refCell) {
+                            float xCell = 0.0;
+                            float yCell = 0.0;
+                            float xEnd = 0.0;
+                            float yEnd = 0.0;
+                            // find coordinate of star on cell
+                            pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y);
+                            // find coordinate of end-point on chip
+
+                            int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE");
+                            pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize);
+
+                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
+                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
+                        }
+                    }
                 }
 
@@ -259 +278 @@
                 // supplied here.
                 psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
-                if (ghosts) { 
-                    // mask the ghosts on this readout
-                    for (int i = 0; i < ghosts->n; i++) {
-                        psastroGhost *ghost = ghosts->data[i];
-                        // XXX bright vs faint ghost bits? (OR with SUSPECT) 
-                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
-                    }
-                }
+                if (ghosts) {
+                    // mask the ghosts on this readout
+                    for (int i = 0; i < ghosts->n; i++) {
+                        psastroGhost *ghost = ghosts->data[i];
+                        // XXX bright vs faint ghost bits? (OR with SUSPECT)
+                        psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
+                    }
+                }
 
                 // Select the glint mask regions for this readout (loaded in
@@ -275 +294 @@
                 psArray *glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS");
                 if (glints) {
-                    // mask the glints on this readout
-                    for (int i = 0; i < glints->n; i++) {
-                        psRegion *glint = glints->data[i];
-                        psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
-                    }
-                }
-
-                // this probably should move into a function of its own:
-                {
-                    // select the raw objects for this readout, flag is they fall in a mask
-                    psArray *inSources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
-                    if (inSources == NULL) continue;
-                
-                    // create a replacement output array:
-                    // psArray *outSources = psAllocArrayEmpty(100);
-
-                    // XXX finish this: raise a bit for stars that land on certain types of masks;
-                    // others (eg, bright star core) should be ignored.
-                    for (int i = 0; i < inSources->n; i++) {
-                        pmSource *source = inSources->data[i];
-
-                        int xChip = source->peak->x;
-                        int yChip = source->peak->y;
-
-                        bool onChip = true;
-                        onChip &= (xChip >= 0);
-                        onChip &= (xChip < readoutMask->mask->numCols);
-                        onChip &= (yChip >= 0);
-                        onChip &= (yChip < readoutMask->mask->numRows);
-                        if (!onChip) {
-                            // if the source is off the edge of the chip, raise a different bit?
-                            source->mode |= PM_SOURCE_MODE_OFF_CHIP;
-                            continue;
-                        }
-
-                        psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
-                        if (value & ghostMaskValue) {
-                            source->mode |= PM_SOURCE_MODE_ON_GHOST;
-                        }
-                        // XXX note that for now, glint and ghost are identical
-                        pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
-                        if (value & glintMaskValue) {
-                            source->mode |= PM_SOURCE_MODE_ON_GLINT;
-                        }
-                        if (value & spikeMaskValue) {
-                            source->mode |= PM_SOURCE_MODE_ON_SPIKE;
-                        }
-                    }
-                }
+                    // mask the glints on this readout
+                    for (int i = 0; i < glints->n; i++) {
+                        psRegion *glint = glints->data[i];
+                        psastroMaskRectangle (readoutMask->mask, glintMaskValue, glint->x0, glint->y0, glint->x1, glint->y1);
+                    }
+                }
+
+                // Select the crosstalk artifact regions for this readout, and mask a circular region
+                // corresponding to the source star
+
+                psArray *crosstalks = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.CROSSTALKS");
+                if (crosstalks) {
+                  for (int i = 0; i < crosstalks->n; i++) {
+                    pmAstromObj *ref = crosstalks->data[i];
+                    float radius = REFSTAR_MASK_CROSSTALK_MAG_SLOPE * (REFSTAR_MASK_CROSSTALK_MAG_MAX - ref->Mag);
+                    psTrace("psastro.crosstalk",2,"Masking star on Chip %s @ (%f,%f) Magnitude: %f Radius %f\n",
+                            psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"),
+                            ref->chip->x,ref->chip->y,ref->Mag,radius);
+                    // XXX for now, assume cell binning is 1x1 relative to chip
+                    psastroMaskCircle (readoutMask->mask, crosstalkMaskValue, ref->chip->x, ref->chip->y, radius, radius);
+                  }
+                }
+
+                // this probably should move into a function of its own:
+                {
+                    // select the raw objects for this readout, flag is they fall in a mask
+                    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+                    if (!detections) continue;
+
+                    psArray *inSources = detections->allSources;
+                    psAssert (inSources, "missing sources?");
+
+                    // create a replacement output array:
+                    // psArray *outSources = psAllocArrayEmpty(100);
+
+                    // XXX finish this: raise a bit for stars that land on certain types of masks;
+                    // others (eg, bright star core) should be ignored.
+                    for (int i = 0; i < inSources->n; i++) {
+                        pmSource *source = inSources->data[i];
+
+                        int xChip = source->peak->x;
+                        int yChip = source->peak->y;
+
+                        bool onChip = true;
+                        onChip &= (xChip >= 0);
+                        onChip &= (xChip < readoutMask->mask->numCols);
+                        onChip &= (yChip >= 0);
+                        onChip &= (yChip < readoutMask->mask->numRows);
+                        if (!onChip) {
+                            // if the source is off the edge of the chip, raise a different bit?
+                            source->mode |= PM_SOURCE_MODE_OFF_CHIP;
+                            continue;
+                        }
+
+                        psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
+                        if (value & ghostMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_GHOST;
+                        }
+                        // XXX note that for now, glint and ghost are identical
+                        pmSourceMode PM_SOURCE_MODE_ON_GLINT = PM_SOURCE_MODE_ON_GHOST;
+                        if (value & glintMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_GLINT;
+                        }
+                        if (value & spikeMaskValue) {
+                            source->mode |= PM_SOURCE_MODE_ON_SPIKE;
+                        }
+                    }
+                }
             }
         }
@@ -344 +382 @@
 
     if (COUNT_GHOSTS) {
-        // save nGhosts to update header.
-        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
-        if (!updates) {
-            updates = psMetadataAlloc ();
-            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
-            psFree (updates);
-        }
-        psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
+        // save nGhosts to update header.
+        psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
+        if (!updates) {
+            updates = psMetadataAlloc ();
+            psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
+            psFree (updates);
+        }
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts);
     }