- Timestamp:
- May 26, 2009, 1:59:32 PM (17 years ago)
- Location:
- branches/cnb_branches/cnb_branch_20090301
- Files:
-
- 4 edited
-
. (modified) (1 prop)
-
psastro (modified) (1 prop)
-
psastro/src (modified) (1 prop)
-
psastro/src/psastroMaskUpdates.c (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/cnb_branches/cnb_branch_20090301
- Property svn:mergeinfo changed
-
branches/cnb_branches/cnb_branch_20090301/psastro
- Property svn:mergeinfo changed
-
branches/cnb_branches/cnb_branch_20090301/psastro/src
- Property svn:ignore
-
old new 17 17 psastroModelFit 18 18 psastroExtract 19 psastroVersionDefinitions.h
-
- Property svn:ignore
-
branches/cnb_branches/cnb_branch_20090301/psastro/src/psastroMaskUpdates.c
r21409 r24244 19 19 } 20 20 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 31 21 /** 32 22 * create a mask or mask regions based on the collection of reference stars that * are in the vicinity of each chip … … 40 30 float zeropt, exptime; 41 31 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 43 35 44 36 // psImageMaskType maskBlank = pmConfigMaskGet("BLANK", config); // Mask value for blank pixels … … 55 47 if (!REFSTAR_MASK) return true; 56 48 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 bool COUNT_GHOSTS = psMetadataLookupF32 (&status, recipe, "REFSTAR_COUNT_GHOSTS"); 55 double GHOST_MAX_MAG = psMetadataLookupF32 (&status, recipe, "GHOST_MAX_MAG"); 56 int nGhosts = 0; 57 57 58 psLogMsg ("psastro", PS_LOG_INFO, "generating a bright-star mask"); 58 59 59 bool REFSTAR_MASK_ REGIONS = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_REGIONS");60 bool REFSTAR_MASK_BLEED = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED"); 60 61 61 62 double REFSTAR_MASK_MAX_MAG = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG"); … … 90 91 REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset; 91 92 REFSTAR_MASK_BLEED_MAG_MAX += MagOffset; 93 GHOST_MAX_MAG += MagOffset; 92 94 93 95 // select the output mask image :: we mosaic to chip mosaic format … … 132 134 133 135 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE; 134 135 // text region files for testing136 FILE *f = NULL;137 if (REFSTAR_MASK_REGIONS) {138 char *filename = NULL;139 char *chipname = psMetadataLookupStr (&status, chip->concepts, "CHIP.NAME");140 psStringAppend (&filename, "refstars.mask.%s.dat", chipname);141 FILE *f = fopen (filename, "w");142 if (!f) {143 psWarning ("cannot create refstar mask file %s\n", filename);144 continue;145 }146 psFree (filename);147 }148 136 149 137 pmChip *refChip = pmFPAviewThisChip (view, refMask->fpa); … … 172 160 173 161 // we mask pixels on the input mask image (chip-mosaic) 174 pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa);175 pmReadout *readoutMask = NULL;176 if (cellMask->readouts->n) {177 readoutMask = cellMask->readouts->data[0];178 }162 // pmCell *cellMask = pmFPAviewThisCell(view, outMask->fpa); 163 // pmReadout *readoutMask = NULL; 164 // if (cellMask->readouts->n) { 165 // readoutMask = cellMask->readouts->data[0]; 166 // } 179 167 180 168 // process each of the readouts … … 183 171 if (! readout->data_exists) { continue; } 184 172 173 // XXX why not do this? 174 pmReadout *readoutMask = pmFPAviewThisReadout (view, outMask->fpa); 175 if (!readoutMask) continue; 176 185 177 // select the raw objects for this readout 186 178 psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); 187 if ( refstars == NULL) { continue; }179 if (!refstars) continue; 188 180 189 181 // we need to generate the following masks regions: … … 191 183 // 2) diffraction spikes in direction ROT - ROTo 192 184 // 3) bleed trail in the direction of the readout 193 194 // XXX for the moment, let's just generate mana region-file objects in chip pixel space 185 // update: with 'burntool' applied to the data, the bleed trail mask is not needed. 186 195 187 for (int i = 0; i < refstars->n; i++) { 196 188 pmAstromObj *ref = refstars->data[i]; 189 190 if (COUNT_GHOSTS) { 191 if (ref->Mag > GHOST_MAX_MAG) { 192 nGhosts ++; 193 } 194 } 195 197 196 if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue; 198 197 199 // XXX convert ref->Mag to instrumental mags 198 // the reference magnitudes have been converted from the instrumental 199 // values supplied in the recipe to apparent mags (above) 200 200 201 201 // CIRCLE around the stars (scaled by magnitude) 202 202 float radius = REFSTAR_MASK_SATSTAR_MAG_SLOPE * (REFSTAR_MASK_SATSTAR_MAG_MAX - ref->Mag); 203 203 204 if (REFSTAR_MASK_REGIONS) { 205 fprintf (f, "CIRCLE %f %f %f %f\n", ref->chip->x, ref->chip->y, radius, radius); 204 // XXX for now, assume cell binning is 1x1 relative to chip 205 psastroMaskCircle (readoutMask->mask, starMaskValue, ref->chip->x, ref->chip->y, radius, radius); 206 207 for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) { 208 209 float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO; 210 211 // LINE for boundaries of the saturation spikes (scaled by magnitude) 212 // float MAG_MAX = (theta == 0.0) || (theta == M_PI) ? REFSTAR_MASK_SATSPIKE_MAG_MAX1 : REFSTAR_MASK_SATSPIKE_MAG_MAX2; 213 214 float MAG_MAX = REFSTAR_MASK_SATSPIKE_MAG_MAX; 215 float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (MAG_MAX - ref->Mag); 216 float spikeWidth = 0.5*REFSTAR_MASK_SATSPIKE_WIDTH; 217 // XXX we can make the width depend on the spike as well... 218 // The length should also be a function of the image background level 219 220 psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); 206 221 } 207 222 208 // XXX for now, assume cell binning is 1x1 relative to chip 209 if (readoutMask) { 210 psastroMaskCircle (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, radius, radius); 211 } 212 213 // LINE for boundaries of the saturation spikes (scaled by magnitude) 214 float spikeLength = REFSTAR_MASK_SATSPIKE_MAG_SLOPE * (REFSTAR_MASK_SATSPIKE_MAG_MAX - ref->Mag); 215 float spikeWidth = REFSTAR_MASK_SATSPIKE_WIDTH; 216 217 for (float theta = 0.0; theta < 2*M_PI; theta += M_PI / 2.0) { 218 float x0, y0, x1, y1, dx, dy; 219 220 float Theta = theta - ROTANGLE - REFSTAR_MASK_SATSTAR_POS_ZERO; 221 222 if (REFSTAR_MASK_REGIONS) { 223 // lower side 224 x0 = ref->chip->x + spikeWidth*sin(Theta); 225 y0 = ref->chip->y - spikeWidth*cos(Theta); 226 x1 = ref->chip->x + spikeLength*cos(Theta) + spikeWidth*sin(Theta); 227 y1 = ref->chip->y + spikeLength*sin(Theta) - spikeWidth*cos(Theta); 228 dx = x1 - x0; 229 dy = y1 - y0; 230 231 fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); 232 233 // upper side 234 x0 = ref->chip->x - spikeWidth*sin(Theta); 235 y0 = ref->chip->y + spikeWidth*cos(Theta); 236 x1 = ref->chip->x + spikeLength*cos(Theta) - spikeWidth*sin(Theta); 237 y1 = ref->chip->y + spikeLength*sin(Theta) + spikeWidth*cos(Theta); 238 dx = x1 - x0; 239 dy = y1 - y0; 240 fprintf (f, "LINE %f %f %f %f\n", x0, y0, dx, dy); 241 } 242 243 if (readoutMask) { 244 psastroMaskBox (readoutMask->mask, maskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); 245 } 246 } 247 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 (REFSTAR_MASK_REGIONS) { 253 fprintf (f, "LINE %f %f %f %f\n", ref->chip->x, ref->chip->y, 0.0, -100.0); 254 } 255 256 if (readoutMask && refCell) { 257 float xCell = 0.0; 258 float yCell = 0.0; 259 float xEnd = 0.0; 260 float yEnd = 0.0; 261 // find coordinate of star on cell 262 pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y); 263 // find coordinate of end-point on chip 264 265 int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE"); 266 pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize); 267 268 float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); 269 psastroMaskRectangle (readoutMask->mask, maskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); 270 } 223 if (REFSTAR_MASK_BLEED) { 224 225 // convert x,y chip coordinates to cells in maskChip 226 pmCell *refCell = pmCellInChip (refChip, ref->chip->x, ref->chip->y); 227 228 // LINE for boundaries of the bleed lines 229 if (refCell) { 230 float xCell = 0.0; 231 float yCell = 0.0; 232 float xEnd = 0.0; 233 float yEnd = 0.0; 234 // find coordinate of star on cell 235 pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y); 236 // find coordinate of end-point on chip 237 238 int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE"); 239 pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize); 240 241 float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); 242 psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); 243 } 244 } 271 245 } 246 247 // select the raw objects for this readout (loaded in psastroExtract.c) 248 psArray *ghosts = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GHOSTS"); 249 if (ghosts == NULL) { continue; } 250 251 // identify the bright stars of interest 252 for (int i = 0; i < ghosts->n; i++) { 253 psastroGhost *ghost = ghosts->data[i]; 254 // XXX bright vs faint ghost bits? (OR with SUSPECT) 255 psastroMaskEllipticalAnnulus (readoutMask->mask, ghostMaskValue, ghost->chip->x, ghost->chip->y, ghost->inner, ghost->outer); 256 } 257 258 // select the raw objects for this readout, flag is they fall in a mask 259 psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS"); 260 if (rawstars == NULL) return false; 261 262 // XXX finish this: 263 for (int i = 0; i < rawstars->n; i++) { 264 pmAstromObj *raw = rawstars->data[i]; 265 psImageMaskType value = readoutMask->mask->data.PS_TYPE_IMAGE_MASK_DATA[(int)(raw->chip->x)][(int)(raw->chip->y)]; 266 if (value) continue; 267 } 272 268 } 273 269 } 274 270 275 // output sequence for mask corresponding to this chip 271 // output sequence for mask corresponding to this chip (XXX this may not be needed...) 276 272 *viewMask = *view; 277 273 while ((cell = pmFPAviewNextCell (viewMask, outMask->fpa, 1)) != NULL) { … … 285 281 if (!pmFPAfileIOChecks (config, viewMask, PM_FPA_AFTER)) ESCAPE; 286 282 } 287 288 283 if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; 289 if (REFSTAR_MASK_REGIONS) {290 fclose (f);291 }292 284 } 293 285 if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE; 286 287 if (COUNT_GHOSTS) { 288 // save nGhosts to update header. 289 psMetadata *updates = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER"); 290 if (!updates) { 291 updates = psMetadataAlloc (); 292 psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", updates); 293 psFree (updates); 294 } 295 psMetadataAddS32 (updates, PS_LIST_TAIL, "NGHOSTS", PS_META_REPLACE, "total expected ghosts", nGhosts); 296 } 294 297 295 298 // deactivate all files … … 301 304 } 302 305 303 // XXX this is going to be very slow...304 pmCell *pmCellInChip (pmChip *chip, float x, float y) {305 306 # if (0)307 int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");308 int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");309 int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");310 int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");311 312 // XXX fix the binning : currently not selected from concepts313 // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y314 // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y315 int xBin = 1;316 int yBin = 1;317 318 // Position on the cell319 float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, xBin);320 float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, yBin);321 # endif322 323 for (int i = 0; i < chip->cells->n; i++) {324 325 pmCell *cell = chip->cells->data[i];326 psRegion *region = pmCellExtent (cell);327 328 if (x < region->x0) goto skip;329 if (x > region->x1) goto skip;330 if (y < region->y0) goto skip;331 if (y > region->y1) goto skip;332 333 psFree (region);334 return cell;335 336 skip:337 psFree (region);338 }339 return NULL;340 }341 342 /**343 * convert chip coords to cell coords, given known cell344 */345 bool pmCellCoordsForChip (float *xCell, float *yCell, pmCell *cell, float xChip, float yChip) {346 347 int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");348 int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");349 int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");350 int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");351 352 // XXX fix the binning : currently not selected from concepts353 // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y354 // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y355 int xBin = 1;356 int yBin = 1;357 358 // Position on the cell359 // ((pos)*(binning)*(cellParity) + (cell0))360 // XXX this is probably totally wrong now....361 // ((pos) - (cell0))*(cellParity)/(binning))362 *xCell = (xChip - x0Cell)*xParityCell/xBin;363 *yCell = (yChip - y0Cell)*yParityCell/yBin;364 365 return true;366 }367 368 /**369 * convert chip coords to cell coords, given known cell370 */371 bool pmChipCoordsForCell (float *xChip, float *yChip, pmCell *cell, float xCell, float yCell) {372 373 int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");374 int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");375 int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");376 int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");377 378 // XXX fix the binning : currently not selected from concepts379 // int xBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN"); // Binning in x and y380 // int yBin = psMetadataLookupS32(NULL, cell->concepts, "CELL.YBIN"); // Binning in x and y381 int xBin = 1;382 int yBin = 1;383 384 // Position on the cell385 // ((pos)*(binning)*(cellParity) + (cell0))386 // XXX this is probably totally wrong now....387 // ((pos) - (cell0))*(cellParity)/(binning))388 *xChip = xCell*xBin*xParityCell + x0Cell;389 *yChip = yCell*yBin*yParityCell + y0Cell;390 391 return true;392 }393 394 // XXX should be doing an OR395 bool psastroMaskCircle (psImage *mask, psImageMaskType value, float x0, float y0, float dX, float dY) {396 397 // XXX need to worry about row0, col0398 for (int ix = -dX; ix <= +dX; ix++) {399 int jx = ix + x0;400 if (jx < 0) continue;401 if (jx >= mask->numCols) continue;402 for (int iy = -dY; iy <= +dY; iy++) {403 int jy = iy + y0;404 if (jy < 0) continue;405 if (jy >= mask->numRows) continue;406 407 double r2 = PS_SQR(ix/dX) + PS_SQR(iy/dY);408 if (r2 > 1.0) continue;409 410 mask->data.PS_TYPE_IMAGE_MASK_DATA[jy][jx] |= value;411 }412 }413 return true;414 }415 416 // XXX should be doing an OR417 bool psastroMaskBox (psImage *mask, psImageMaskType value, float x0, float y0, float dL, float dW, float theta) {418 419 // draw a series of lines (from -0.5*dW to +0.5*dW) of length dL, starting at x0, y0, angle theta420 421 float xs = x0;422 float ys = y0;423 424 float xe = xs + dL*cos(theta);425 float ye = ys + dL*sin(theta);426 427 psastroMaskLine (mask, value, xs, ys, xe, ye, (int) dW);428 return true;429 }430 431 /**432 * identify the quadrant and draw the correct line433 */434 void psastroMaskLine (psImage *mask, psImageMaskType value, double x1, double y1, double x2, double y2, int dW) {435 436 int FlipDirect, FlipCoords;437 int X1, Y1, X2, Y2, dX, dY;438 439 /* rather than draw the line from float positions, we find the closest440 integer end-points and draw the line between those pixels */441 442 X1 = ROUND(x1);443 Y1 = ROUND(y1);444 X2 = ROUND(x2);445 Y2 = ROUND(y2);446 447 dX = X2 - X1;448 dY = Y2 - Y1;449 450 FlipCoords = (abs(dX) < abs(dY));451 FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2);452 453 if (!FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X1, Y1, X2, Y2, dW, FALSE);454 if ( FlipDirect && !FlipCoords) psastroMaskLineBresen (mask, value, X2, Y2, X1, Y1, dW, FALSE);455 if (!FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y1, X1, Y2, X2, dW, TRUE);456 if ( FlipDirect && FlipCoords) psastroMaskLineBresen (mask, value, Y2, X2, Y1, X1, dW, TRUE);457 458 return;459 }460 461 /**462 * use the Bresenham line drawing technique463 * integer-only Bresenham line-draw version which is fast464 */465 void psastroMaskLineBresen (psImage *mask, psImageMaskType value, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) {466 467 int X, Y, dX, dY;468 int e, e2;469 470 dX = X2 - X1;471 dY = Y2 - Y1;472 473 Y = Y1;474 e = 0;475 for (X = X1; X <= X2; X++) {476 if (X > 0) {477 if (swapcoords) {478 if (X >= mask->numRows) continue;479 for (int y = Y - dW; y <= Y + dW; y++) {480 if (y < 0) continue;481 if (y >= mask->numCols) continue;482 mask->data.PS_TYPE_IMAGE_MASK_DATA[X][y] |= value;483 }484 } else {485 if (X >= mask->numCols) continue;486 for (int y = Y - dW; y <= Y + dW; y++) {487 if (y < 0) continue;488 if (y >= mask->numRows) continue;489 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][X] |= value;490 }491 }492 }493 e += dY;494 e2 = 2 * e;495 if (e2 > dX) {496 Y++;497 e -= dX;498 }499 if (e2 < -dX) {500 Y--;501 e += dX;502 }503 }504 return;505 }506 507 void psastroMaskRectangle (psImage *mask, psImageMaskType value, int x0, int y0, int x1, int y1) {508 509 int xs = PS_MAX (0, PS_MIN (mask->numCols, PS_MIN (x0, x1)));510 int xe = PS_MAX (0, PS_MIN (mask->numCols, PS_MAX (x0, x1)));511 int ys = PS_MAX (0, PS_MIN (mask->numRows, PS_MIN (y0, y1)));512 int ye = PS_MAX (0, PS_MIN (mask->numRows, PS_MAX (y0, y1)));513 514 for (int iy = ys; iy < ye; iy++) {515 for (int ix = xs; ix < xe; ix++) {516 mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= value;517 }518 }519 }520
Note:
See TracChangeset
for help on using the changeset viewer.
