IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 10, 2009, 3:00:23 PM (17 years ago)
Author:
eugene
Message:

add masking for ghosts based on model (polynomials in a metadata config file); move mask support functions to psastroMaskUtils; plug leak in psastro

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroMaskUpdates.c

    r23665 r23810  
    1919}
    2020
    21 pmCell *pmCellInChip (pmChip *chip, float x, float y);
    22 bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip);
    23 bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell);
    24 
    25 bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY);
    26 bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta);
    27 void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW);
    28 void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
    29 void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1);
    30 
    3121/**
    3222 * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip
    … …  
    4030    float zeropt, exptime;
    4131
    42     psImageMaskType maskValue  = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
     32    psImageMaskType ghostMaskValue = pmConfigMaskGet("GHOST", config); // Mask value for ghost pixels
     33    psImageMaskType spikeMaskValue = pmConfigMaskGet("SPIKE", config); // Mask value for ghost pixels
     34    psImageMaskType starMaskValue  = pmConfigMaskGet("STARCORE", config); // Mask value for ghost pixels
    4335
    4436    // psImageMaskType maskBlank  = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels
    … …  
    5547    if (!REFSTAR_MASK) return true;
    5648
     49    // convert star positions to ghost positions and add to the readout->analysis data
     50    if (!psastroLoadGhosts (config)) {
     51        psError(PSASTRO_ERR_CONFIG, false, "Error loading ghosts");
     52        return false;
     53    }
     54
    5755    psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask");
    5856
    59     bool REFSTAR_MASK_REGIONS              = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_REGIONS");
    6057    bool REFSTAR_MASK_BLEED                = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED");
    6158
    … …  
    134131        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE;
    135132
    136         // text region files for testing
    137         FILE *f = NULL;
    138         if (REFSTAR_MASK_REGIONS) {
    139           char *filename = NULL;
    140           char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");
    141           psStringAppend (&filename, "refstars.mask.%s.dat", chipname);
    142           FILE *f = fopen (filename, "w");
    143           if (!f) {
    144             psWarning ("cannot create refstar mask file %s\n", filename);
    145             continue;
    146           }
    147           psFree (filename);
    148         }
    149 
    150133        pmChip *refChip  = pmFPAviewThisChip (view, refMask->fpa);
    151134
    … …  
    173156
    174157            // we mask pixels on the input mask image (chip-mosaic)
    175             pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);
    176             pmReadout *readoutMask = NULL;
    177             if (cellMask->readouts->n) {
    178                 readoutMask = cellMask->readouts->data[0];
    179             }
     158            // pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);
     159            // pmReadout *readoutMask = NULL;
     160            // if (cellMask->readouts->n) {
     161            //     readoutMask = cellMask->readouts->data[0];
     162            // }
    180163
    181164            // process each of the readouts
    … …  
    184167                if (! readout->data_exists) { continue; }
    185168
     169                // XXX why not do this?
     170                pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa);
     171                if (!readoutMask) continue;
     172
    186173                // select the raw objects for this readout
    187174                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    188                 if (refstars == NULL) { continue; }
     175                if (!refstars) continue;
    189176
    190177                // we need to generate the following masks regions:
    … …  
    204191                    float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag);
    205192
    206                     if (REFSTAR_MASK_REGIONS) {
    207                       fprintf (f, "CIRCLE %f %f  %f %f\n", ref->chip->x, ref->chip->y, radius, radius);
    208                     }
    209 
    210193                    // XXX for now, assume cell binning is 1x1 relative to chip
    211                     if (readoutMask) {
    212                         psastroMaskCircle (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, radius, radius);
    213                     }
    214 
    215                     // LINE for boundaries of the saturation spikes (scaled by magnitude)
    216                     float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (REFSTAR_MASK_SATSPIKE_MAG_MAX - ref->Mag);
    217                     float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
     194                    psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius);
    218195
    219196                    for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) {
    220                         float x0, y0, x1, y1, dx, dy;
    221197
    222198                        float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO;
    223199
    224                         if (REFSTAR_MASK_REGIONS) {
    225                           // lower side
    226                           x0 = ref->chip->x + spikeWidth*sin(Theta);
    227                           y0 = ref->chip->y - spikeWidth*cos(Theta);
    228                           x1 = ref->chip->x + spikeLength*cos(Theta) + spikeWidth*sin(Theta);
    229                           y1 = ref->chip->y + spikeLength*sin(Theta) - spikeWidth*cos(Theta);
    230                           dx = x1 - x0;
    231                           dy = y1 - y0;
    232 
    233                           fprintf (f, "LINE %f %f  %f %f\n", x0, y0, dx, dy);
    234 
    235                           // upper side
    236                           x0 = ref->chip->x - spikeWidth*sin(Theta);
    237                           y0 = ref->chip->y + spikeWidth*cos(Theta);
    238                           x1 = ref->chip->x + spikeLength*cos(Theta) - spikeWidth*sin(Theta);
    239                           y1 = ref->chip->y + spikeLength*sin(Theta) + spikeWidth*cos(Theta);
    240                           dx = x1 - x0;
    241                           dy = y1 - y0;
    242                           fprintf (f, "LINE %f %f  %f %f\n", x0, y0, dx, dy);
    243                         }
    244 
    245                         if (readoutMask) {
    246                             psastroMaskBox (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
    247                         }
     200                        // LINE for boundaries of the saturation spikes (scaled by magnitude)
     201                        // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2;
     202
     203                        float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX;
     204                        float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag);
     205                        float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH;
     206                        // XXX we can make the width depend on the spike as well...
     207                        // The length should also be a function of the image background level
     208
     209                        psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta);
    248210                    }
    249211
    … …  
    254216
    255217                        // LINE for boundaries of the bleed lines
    256                         if (REFSTAR_MASK_REGIONS) {
    257                             fprintf (f, "LINE %f %f  %f %f\n", ref->chip->x, ref->chip->y, 0.0, -100.0);
    258                         }
    259 
    260                         if (readoutMask && refCell) {
     218                        if (refCell) {
    261219                            float xCell = 0.0;
    262220                            float yCell = 0.0;
    … …  
    271229
    272230                            float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag);
    273                             psastroMaskRectangle (readoutMask->mask, maskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
     231                            psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd);
    274232                        }
    275233                    }
    276234                }
     235
     236                // select the raw objects for this readout (loaded in psastroExtract.c)
     237                psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS");
     238                if (ghosts == NULL) { continue; }
     239
     240                // identify the bright stars of interest
     241                for (int i = 0; i < ghosts->n; i++) {
     242                    psastroGhost *ghost = ghosts->data[i];
     243                    // XXX bright vs faint ghost bits? (OR with SUSPECT)
     244                    psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer);
     245                }
     246
     247                // select the raw objects for this readout, flag is they fall in a mask
     248                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
     249                if (rawstars == NULL) return false;
     250               
     251                // XXX finish this:
     252                for (int i = 0; i < rawstars->n; i++) {
     253                    pmAstromObj *raw = rawstars->data[i];
     254                    psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)];
     255                    if (value) continue;
     256                }
    277257            }
    278258        }
    279259
    280         // output sequence for mask corresponding to this chip
     260        // output sequence for mask corresponding to this chip (XXX this may not be needed...)
    281261        *viewMask = *view;
    282262        while ((cell = pmFPAviewNextCell (viewMask, outMask->fpa, 1)) != NULL) {
    … …  
    290270            if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE;
    291271        }
    292 
    293272        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    294         if (REFSTAR_MASK_REGIONS) {
    295           fclose (f);
    296         }
    297273    }
    298274    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE;
    … …  
    306282}
    307283
    308 // XXX this is going to be very slow...
    309 pmCell *pmCellInChip (pmChip *chip, float x, float y) {
    310 
    311 # if (0)
    312     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    313     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    314     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    315     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    316 
    317     // XXX fix the binning : currently not selected from concepts
    318     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    319     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    320     int xBin = 1;
    321     int yBin = 1;
    322 
    323     // Position on the cell
    324     float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);
    325     float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);
    326 # endif
    327 
    328     for (int i = 0; i < chip->cells->n; i++) {
    329 
    330         pmCell *cell = chip->cells->data[i];
    331         psRegion *region = pmCellExtent (cell);
    332 
    333         if (x < region->x0) goto skip;
    334         if (x > region->x1) goto skip;
    335         if (y < region->y0) goto skip;
    336         if (y > region->y1) goto skip;
    337 
    338         psFree (region);
    339         return cell;
    340 
    341     skip:
    342         psFree (region);
    343     }
    344     return NULL;
    345 }
    346 
    347 /**
    348  * convert chip coords to cell coords, given known cell
    349  */
    350 bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip) {
    351 
    352     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    353     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    354     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    355     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    356 
    357     // XXX fix the binning : currently not selected from concepts
    358     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    359     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    360     int xBin = 1;
    361     int yBin = 1;
    362 
    363     // Position on the cell
    364     // ((pos)*(binning)*(cellParity) + (cell0))
    365     // XXX this is probably totally wrong now....
    366     // ((pos) - (cell0))*(cellParity)/(binning))
    367     *xCell = (xChip - x0Cell)*xParityCell/xBin;
    368     *yCell = (yChip - y0Cell)*yParityCell/yBin;
    369 
    370     return true;
    371 }
    372 
    373 /**
    374  * convert chip coords to cell coords, given known cell
    375  */
    376 bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell) {
    377 
    378     int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    379     int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    380     int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    381     int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    382 
    383     // XXX fix the binning : currently not selected from concepts
    384     // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y
    385     // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y
    386     int xBin = 1;
    387     int yBin = 1;
    388 
    389     // Position on the cell
    390     // ((pos)*(binning)*(cellParity) + (cell0))
    391     // XXX this is probably totally wrong now....
    392     // ((pos) - (cell0))*(cellParity)/(binning))
    393     *xChip = xCell*xBin*xParityCell + x0Cell;
    394     *yChip = yCell*yBin*yParityCell + y0Cell;
    395 
    396     return true;
    397 }
    398 
    399 // XXX should be doing an OR
    400 bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY) {
    401 
    402     // XXX need to worry about row0, col0
    403     for (int ix = -dX; ix <= +dX; ix++) {
    404         int jx = ix + x0;
    405         if (jx < 0) continue;
    406         if (jx >= mask->numCols) continue;
    407         for (int iy = -dY; iy <= +dY; iy++) {
    408             int jy = iy + y0;
    409             if (jy < 0) continue;
    410             if (jy >= mask->numRows) continue;
    411 
    412             double r2 = PS_SQR(ix/dX) + PS_SQR(iy/dY);
    413             if (r2 > 1.0) continue;
    414 
    415             mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value;
    416         }
    417     }
    418     return true;
    419 }
    420 
    421 // XXX should be doing an OR
    422 bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta) {
    423 
    424     // draw a series of lines (from -0.5*dW to +0.5*dW) of length dL, starting at x0, y0, angle theta
    425 
    426     float xs = x0;
    427     float ys = y0;
    428 
    429     float xe = xs + dL*cos(theta);
    430     float ye = ys + dL*sin(theta);
    431 
    432     psastroMaskLine (mask, value, xs, ys, xe, ye, (int) dW);
    433     return true;
    434 }
    435 
    436 /**
    437  * identify the quadrant and draw the correct line
    438  */
    439 void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW) {
    440 
    441   int FlipDirect, FlipCoords;
    442   int X1, Y1, X2, Y2, dX, dY;
    443 
    444   /* rather than draw the line from float positions, we find the closest
    445      integer end-points and draw the line between those pixels */
    446 
    447   X1 = ROUND(x1);
    448   Y1 = ROUND(y1);
    449   X2 = ROUND(x2);
    450   Y2 = ROUND(y2);
    451 
    452   dX = X2 - X1;
    453   dY = Y2 - Y1;
    454 
    455   FlipCoords = (abs(dX) < abs(dY));
    456   FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2);
    457 
    458   if (!FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X1, Y1, X2, Y2, dW, FALSE);
    459   if ( FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X2, Y2, X1, Y1, dW, FALSE);
    460   if (!FlipDirect &&  FlipCoords) psastroMaskLineBresen (mask, value, Y1, X1, Y2, X2, dW, TRUE);
    461   if ( FlipDirect &&  FlipCoords) psastroMaskLineBresen (mask, value, Y2, X2, Y1, X1, dW, TRUE);
    462 
    463   return;
    464 }
    465 
    466 /**
    467  * use the Bresenham line drawing technique
    468  * integer-only Bresenham line-draw version which is fast
    469  */
    470 void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) {
    471 
    472     int X, Y, dX, dY;
    473     int e, e2;
    474 
    475     dX = X2 - X1;
    476     dY = Y2 - Y1;
    477 
    478     Y = Y1;
    479     e = 0;
    480     for (X = X1; X <= X2; X++) {
    481         if (X > 0) {
    482             if (swapcoords) {
    483                 if (X >= mask->numRows) continue;
    484                 for (int y = Y - dW; y <= Y + dW; y++) {
    485                     if (y < 0) continue;
    486                     if (y >= mask->numCols) continue;
    487                     mask->data.PS_TYPE_IMAGE_MASK_DATA[X][y] |= value;
    488                 }
    489             } else {
    490                 if (X >= mask->numCols) continue;
    491                 for (int y = Y - dW; y <= Y + dW; y++) {
    492                     if (y < 0) continue;
    493                     if (y >= mask->numRows) continue;
    494                     mask->data.PS_TYPE_IMAGE_MASK_DATA[y][X] |= value;
    495                 }
    496             }
    497         }
    498         e += dY;
    499         e2 = 2 * e;
    500         if (e2 > dX) {
    501             Y++;
    502             e -= dX;
    503         }
    504         if (e2 < -dX) {
    505             Y--;
    506             e += dX;
    507         }
    508     }
    509     return;
    510 }
    511 
    512 void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1) {
    513 
    514     int xs = PS_MAX (0, PS_MIN (mask->numCols, PS_MIN (x0, x1)));
    515     int xe = PS_MAX (0, PS_MIN (mask->numCols, PS_MAX (x0, x1)));
    516     int ys = PS_MAX (0, PS_MIN (mask->numRows, PS_MIN (y0, y1)));
    517     int ye = PS_MAX (0, PS_MIN (mask->numRows, PS_MAX (y0, y1)));
    518 
    519     for (int iy = ys; iy < ye; iy++) {
    520         for (int ix = xs; ix < xe; ix++) {
    521             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= value;
    522         }
    523     }
    524 }
    525 
Note: See TracChangeset for help on using the changeset viewer.