Changeset 39232 for trunk/doc/release.2015/ps1.detrend/detrend.tex
- Timestamp:
- Dec 4, 2015, 5:33:06 PM (11 years ago)
- File:
-
- 1 edited
-
trunk/doc/release.2015/ps1.detrend/detrend.tex (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/doc/release.2015/ps1.detrend/detrend.tex
r38557 r39232 139 139 140 140 \section{INTRODUCTION}\label{sec:intro} 141 141 %% http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?2007ASPC..364..153M&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf 142 142 \section{Camera description} 143 143 144 \czwdraft{reference to original paper} 145 146 \czwdraft{60 otas} 147 148 \czwdraft{64 cells per ota} 149 150 \czwdraft{effectively 60x64 different cameras, each with particular gain/noise/etc characteristics} 151 152 \czwdraft{Add summary of detrending steps} 153 154 \czwdraft{Summary of detrending steps with references to the sections} 144 The Pan-STARRS 1 Science Survey uses the 1.4 giga-pixel GPC1 camera with the PS1 telescope on Haleakala Maui to image the sky north of $-30$ declination. The GPC1 camera is composed of 60 orthogonal transfer array (OTA) devices, each of with is an $8\times{}8$ grid of readout cells. This parallelizes the readout process, reducing the overhead in each exposure. However, as a consequence of this large number of individual detector readouts, there are a number of calibrations that need to be included to ensure the response is the same across the entire field of view. 145 146 The Pan-STARRS image processing pipeline (IPP) is described elsewhere \citep{MagnierXXX}, but a short summary follows. The archive of raw exposures is stored on disk, with a database storing the metadata of exposure parameters. For the PV3 processing, large contiguous regions were defined, and the images for all exposures within that region lauched for the CHIP stage processing. This stage performs the image detrending (described below in section \ref{dead ref}), as well as the single epoch photometry \citep{MagnierXXY}. Following the CHIP stage is the CAMERA stage, in which the astrometry and photometry for entire exposure is calibrated against the reference catalog. This stage also performs masking updates based on the now-known positions and brightnesses of stars that create dynamic features (see \ref{dynamic_masks} below). The WARP stage is the next to operate on the data, transforming the detector oriented CHIP stage images into sky-oriented images that have common tesselations and sky projections (Section \ref{warping}). When all WARP stage processing is done for a region of the sky, STACK processing is performed (Section \ref{stacking}) to construct deeper, fully populated images from the set of WARP images that cover that region of the sky. Beyond the STACK stage, a series of addition stages are done that are described in other papers. Transient features are identified in the DIFF stage, which takes input WARP and/or STACK data and performs image differencing \citep{HuberXXX}. Further photometry is performed in the STATICSKY and SKYCAL stages, which add extended source fitting to the point source photometry of objects detected in the STACK images, and calibrate the results against the reference catalog. The FULLFORCE stage takes the catalog output of the SKYCAL stage, and uses the objects detected in that to perform forced photometry on the individual WARP stage images. The details of this photometry are provided in \citet{MagnierXXY}. 147 148 The full detrend application and processing are described in detail in the the sections below, but a short summary follows. Once an exposure has been observed on the summit, it is transferred to the main IPP processing cluster at the MRTC-B and registered into the processing database. This triggers a new chip stage reduction for each of the 60 OTA images that detrends and mosaicks the individual readout cells before measuring the photometric properties of the astronomical objects detected therein. To begin the detrending, pre-determined static bad pixel masks are used to exclude detector regions that are known to be uncorrectable. Following this, persisitence trails related to the incomplete transfer of charge in the readout process are corrected. The image overscan is subtracted, and the known non-linearity of full readouts cells and the substantially worse issues on the edges of the cells are boosted to the expected levels. The temperature and exposure time dependent dark model is then applied, and the noisemap related to the correlated read noise that is not fully corrected by the dark model is then calculated. The flat field correction is applied next, and any fringe correction necessary for long wavelength data is subtracted. Finally, GPC1 specific ``pattern'' corrections are applied to attempt to reduce the cell to cell differences within a single OTA. 149 150 151 152 % Discuss 2-phase/3-phase device differnces 155 153 156 154 \section{Burntool / Persistence effect} 157 155 158 Stars that are nearing saturation \czwdraft{(30000 DN)}cause156 Stars that are nearing saturation on GPC1 cause 159 157 persistance problems during the read out of the image, creating trails 160 158 of light are left on the image. During the read out process of an … … 179 177 BURNTOOL program. This program does an initial scan of the images, 180 178 and identifies stars brighter than a given threshold of 30000 DN. The 181 trail from that star is fit with a one-dimensional power law 182 \czwdraft{in each pixel column}, based on empirical evidence that this 179 trail from that star is fit with a one-dimensional power law in each pixel column, based on empirical evidence that this 183 180 is the functional form of this persistence effect. Once this fit is 184 181 done, the model is subtracted from the image, and the location of the … … 218 215 219 216 \section{Masking} 217 \czwdraft{Technically, we mask the image prior to burntool application now.} 220 218 221 219 \subsection{Static Masks} 222 220 223 221 Due to the large size of the detector, it is to be expected that there 224 will be a number of pixel defects that \czwdraft{do not measure light}225 as well astheir neighbors. To remove these pixels, we have222 will be a number of pixel defects that do not have the detection sensitivity on par 223 with their neighbors. To remove these pixels, we have 226 224 constructed a static mask that identifies the known defects. This 227 225 mask is constructed in three phases. … … 231 229 detector. Twenty-five of the sixty OTAs in GPC1 show some evidence of 232 230 CTE issues, with this pattern showing up (to varying degrees) in 233 triangular sets of cells on the OTA. \czwdraft{probably a figure would 234 help explain this?} To generate the mask, a sample set of evenly 235 illuminated flat field images are measured to produce a map of the 236 image variance in 20x20 pixel bins. As the flat image largely 237 illuminates the image uniformly, the expected variances should be 238 Poissonian distributed with the flux level. However, in regions with 239 CTE issues, adjacent pixels are not independent, allowing the charge 240 to spread. This reduces the pixel-to-pixel differences, resulting in 241 a lower-than-expected variance. All regions with variance 242 \czwdraft{X} smaller than expected are added to the static CTEMASK. 243 244 The next step of mask construction is to examine the detector for 231 roughly triangular patches on the OTA due to defects in the 232 semiconductor \czwdraft{doping}. To generate the mask for these 233 regions, a sample set of \czwdraft{N} evenly illuminated flat field 234 images were measured to produce a map of the image variance in 20x20 235 pixel bins. As the flat image is expected to illuminate the image 236 uniformly, the expected variances in each bin should be Poissonian 237 distributed with the flux level. However, in regions with CTE issues, 238 adjacent pixels are not independent, allowing the charge in those 239 pixels to spread. This reduces the pixel-to-pixel differences, 240 resulting in a lower-than-expected variance. All regions with 241 variance \czwdraft{0.5} smaller than expected are added to the static 242 CTEMASK. 243 244 The next step of mask construction is to examine the flat and dark 245 models, and exclude pixels that appear to be poorly corrected by these 246 models. The darkmask process looks for pixels that are more than 247 \czwdraft{8} sigma discrepant in \czwdraft{10\%} of the 248 \czwdraft{test} images after those images have had the dark model 249 applied to them. These pixels are assumed to be unstable with respect 250 to the dark model, and have the DARK bit set in the static mask, 251 indicating that they are unreliable in scientific observing. 252 Similarly, the flatmask process looks for pixels that are \czwdraft{3} 253 sigma discrepant in the same fraction of \czwdraft{test} images after 254 both the dark and flat models have been applied. Those pixels that do 255 not follow the flat field model of the rest of image are assigned the 256 FLAT mask bit in the static mask, removing the pixels that cannot be 257 corrected to a linear response. 258 259 The final step of mask construction is to examine the detector for 245 260 bright columns and other static pixel issues. This is first done by 246 \czwdraft{I think Heather wrote a program to do this, but I'm not 247 totally sure how it works} scanning a set of images for pixels that 248 have values that do not change throughout a sequence of \czwdraft{N} 249 exposures. Such common pixel values cannot be caused by astronomical 250 effects, and must be due to the detector itself. This does an 251 excellent job of removing the majority of the problem pixels. A 252 manual inspection allows human interaction to identify other 253 inconsistent pixels including the vignetted regions around the edge of 254 the detector. \czwdraft{This might be a lie} As the size of the 255 vignetted region changes with filter, we have taken the g filter as 256 the baseline to define the static mask, resulting in the smallest 257 possible unvignetted region. 258 259 The final static mask is the union of the CTE mask, the manual mask, \czwdraft{make this a paragraph}. 261 processing a set of \czwdraft{100 i filter} science images in the same 262 fashion as for the darktest. A median image is constructed from these 263 inputs along with the per-pixel variance. These images are used to 264 identify pixels that have unexpectedly low variation between all 265 inputs, as well as those that significantly deviate from the global 266 median value. Once this initial set of bad pixels is identified, a 267 $3\times{}3$ pixel triangular kernel is convolved with the initial 268 set, and any convolved pixel with value greater than \czwdraft{1.0} is 269 assigned to the static mask. This does an excellent job of removing 270 the majority of the problem pixels. A subsequent manual inspection 271 allows human interaction to identify other inconsistent pixels 272 including the vignetted regions around the edge of the detector. 273 \czwdraft{This might be a lie} As the size of the vignetted region 274 changes with filter, we have used the g filter to set the baseline 275 unvignetted region to define the static mask, resulting in the 276 smallest possible unvignetted region. 260 277 261 278 \begin{figure} … … 268 285 defects, we also generate a set of dynamic masks that change with the 269 286 astronomical features in the image. These masks are advisory in 270 nature, and no not completely exclude the pixel from further 271 consideration. The first of these dynamic masks indicates the 272 presence of a corrected burntool trail. These pixels are included for 273 phtometry, but are rejected more readily in the stacking and 274 difference image construction. 287 nature, and do not completely exclude the pixel from further 288 processing consideration. The first of these dynamic masks indicates 289 the presence of a corrected burntool trail. These pixels are included 290 for phtometry, but are rejected more readily in the stacking and 291 difference image construction, as they are more likely to have small 292 residual contributions from the under or over subtraction of the 293 burntool correction. 275 294 276 295 The remaining dynamic masks are not generated until the IPP camera … … 278 297 photometry is complete, and an astrometric solution is known for the 279 298 exposure. This added information provides the positions of bright 280 sources, which are the origin for the image artifacts that the dynamic 281 mask identifies. 299 sources based on the reference catalog, including those that fall 300 slightly out of the detector field of view or within the inter chip 301 gaps, where internal photometry may not have identified them. These 302 bright sources are the origin for many of the image artifacts that the 303 dynamic mask identifies and excludes. 304 305 306 \begin{deluxetable}{ccl} 307 \tablecolumns{3} 308 \tablewidth{0pc} 309 \tablecaption{GPC1 Mask Values} 310 \tablehead{\colhead{Mask Name} & \colhead{Mask Value} & \colhead{Description}} 311 \startdata 312 DETECTOR & 0x0001 & A detector defect is present. \\ 313 FLAT & 0x0002 & The flat field model does not calibrate the pixel reliably. \\ 314 DARK & 0x0004 & The dark model does not calibrate the pixel reliably. \\ 315 BLANK & 0x0008 & The pixel does not contain valid data. \\ 316 CTE & 0x0010 & The pixel has poor charge transfer efficiency. \\ 317 SAT & 0x0020 & The pixel is saturated. \\ 318 LOW & 0x0040 & The pixel has a lower value than expected. \\ 319 SUSPECT & 0x0080 & The pixel is suspected of being bad. \\ 320 BURNTOOL & 0x0080 & The pixel may contain an uncorrected or over-corrected burntool streak. \\ 321 CR & 0x0100 & A cosmic ray is present. \\ 322 SPIKE & 0x0200 & A diffraction spike is present. \\ 323 GHOST & 0x0400 & An optical ghost is present. \\ 324 STREAK & 0x0800 & A streak is present. \\ 325 STARCORE & 0x1000 & A bright star core is present. \\ 326 CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\ 327 CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\ 328 MARK & 0x8000 & An internal flag for temporarily marking a pixel. \\ 329 \enddata 330 \label{tab:mask_values} 331 \end{deluxetable} 332 282 333 283 334 \subsubsection{Crosstalk ghosts} 335 284 336 Due to electrical crosstalk between the flex cables connecting the 285 individual detectors, ghost objects can be created on some OTAs due to 286 the presence of a bright object in a different position. Table 287 \ref{tab:crosstalk_rules} summarizes the list of known crosstalk 288 rules. In each of these cases, a source object brighter than -14.47 289 magnitude (instrumental) creates a ghost object many orders of 290 magnitude fainter at the target location. The cell (x,y) coordinate 291 is identical between source and ghost, as a result of the transfer 292 occurring as the devices are read. A circular mask is asdded to the 293 ghost location with radius $R = 3.44 \left(-14.47 - m_{source, 294 instrumental}\right)$. Any objects in the photometric catalog found 295 at the location of the ghost mask have a flag set, marking the object 296 as a ghost. 297 298 \draft{We also have to deal with bleed ghosts. MAG_MAX = -15, SLOPE = 5.0. Main CT rules only. same OTA Xt=X,Yt=Y,Vt=V,Ut=0:8. width = 5 * (-15 - mag). Top to bottom.} 337 individual detector devices, ghost objects can be created on some OTAs 338 due to the presence of a bright source at a different position on the 339 camera. Table \ref{tab:crosstalk_rules} summarizes the list of known 340 crosstalk rules. In each of these cases, a source object brighter 341 than -14.47 magnitude (instrumental) creates a ghost object many 342 orders of magnitude fainter at the target location. The cell (x,y) 343 coordinate is identical between source and ghost, as a result of the 344 transfer occurring as the devices are read. A circular mask is asdded 345 to the ghost location with radius $R = 3.44 \left(-14.47 - m_{source, 346 instrumental}\right)$ pixels. Any objects in the photometric 347 catalog found at the location of the ghost mask have a \czwdraft{flag} 348 set, marking the object as a likely ghost. The majority of the 349 crosstalk rules are bi-directional, with a source in either position 350 creating a ghost at the corresponding crosstalk target position. The 351 two faintest rules are uni-directional, likely due to differences in 352 the \czwdraft{magical properties of the electronics}. 353 354 For the very brightest sources ($m_{instrumental} < -15$), there can 355 be crosstalk ghosts between all columns of cells during the readout. 356 These ``bleed'' ghosts were originally identified as ghosts of the 357 saturation bleeds appearing in the neighboring cells, and as such, the 358 masking for these objects puts a rectangular mask down from top to 359 bottom of cells in all columns that are in the same row of cells as 360 the bright source. The width of this box is a function of the source 361 magnitude, with $W = 5 * \left(-15 - m_{source, instrumental}\right)$ 362 pixels. 299 363 300 364 \begin{deluxetable}{lllc} … … 318 382 \end{deluxetable} 319 383 384 \begin{figure} 385 \caption{Figure of crosstalk ghost and bright star source. Plot of cut across ghost to illustrate the flat-top shape.} 386 \end{figure} 320 387 321 388 \subsubsection{Optical ghosts} 322 323 Due to an issue with the anti-reflective coating, bright sources can324 also result in large out of focus objects, particularly in the389 % http://arxiv.org/pdf/1207.2513v1.pdf 390 Due to imperfections in the anti-reflective coating, bright sources 391 can also result in large out of focus objects, particularly in the 325 392 g-filter data. These objects are the result of light reflecting back 326 off the surface of the detector, reflecting again off the \czwdraft{No 327 clue} mirror, and then back down onto the focal plane. Due to the 328 extra travel distance, the resulting source is out of focus and 329 elongated along the radial direction of the telescope. These optical 330 ghosts can be modeled as a bright star in location (X,Y) on the focal 331 plane creates a reflection ghost on the opposite side of the optical 332 axis at (-X,-Y). The exact location is fit as a third order 333 polynomial in the focal plane x and y directions. An elliptical 334 annulus mask is constructed at the expected ghost location, with the 335 major and minor axes defined by linear functions of the ghost distance 336 from the optical axis, and orientation \czwdraft{pointing along 337 radius}. All stars brighter than a filter-dependent threshold 338 (listed in table \ref{tab:ghost_magnitudes}) have masks constructed. 393 off the surface of the detector, reflecting again off the lower 394 surfaces of the optics (particularly the L1 corrector lens), and then 395 back down onto the focal plane. Due to the extra travel distance, the 396 resulting source is out of focus and elongated along the radial 397 direction of the telescope. These optical ghosts can be modeled as a 398 bright star in location (X,Y) on the focal plane creates a reflection 399 ghost on the opposite side of the optical axis at (-X,-Y). The exact 400 location is fit as a third order polynomial in the focal plane x and y 401 directions. An elliptical annulus mask is constructed at the expected 402 ghost location, with the major and minor axes defined by linear 403 functions of the ghost distance from the optical axis, and oriented 404 along the radius of the detector. All stars brighter than a 405 filter-dependent threshold (listed in table 406 \ref{tab:ghost_magnitudes}) have such masks constructed. 339 407 340 408 \begin{deluxetable}{lc} … … 356 424 \czwdraft{include full polynomial forms? How best to do that?} 357 425 358 359 %% 360 %% GHOST.CENTER.X METADATA 361 %% NORDER_X S32 3 362 %% NORDER_Y S32 3 363 %% VAL_X00_Y00 F64 -1.215661e+02 364 %% VAL_X01_Y00 F64 1.321875e-02 365 %% VAL_X02_Y00 F64 -4.017026e-09 366 %% VAL_X03_Y00 F64 1.148288e-10 367 %% VAL_X00_Y01 F64 -1.908074e-03 368 %% VAL_X01_Y01 F64 8.479150e-08 369 %% VAL_X02_Y01 F64 1.635732e-11 370 %% VAL_X00_Y02 F64 2.625405e-08 371 %% VAL_X01_Y02 F64 1.125586e-10 372 %% VAL_X00_Y03 F64 2.912432e-12 373 %% NELEMENTS S32 10 374 %% END 375 376 %% GHOST.CENTER.Y METADATA 377 %% NORDER_X S32 3 378 %% NORDER_Y S32 3 379 %% VAL_X00_Y00 F64 2.422174e+01 380 %% VAL_X01_Y00 F64 4.170486e-04 381 %% VAL_X02_Y00 F64 -1.934260e-08 382 %% VAL_X03_Y00 F64 -1.173657e-12 383 %% VAL_X00_Y01 F64 1.189352e-02 384 %% VAL_X01_Y01 F64 -9.256748e-08 385 %% VAL_X02_Y01 F64 1.140772e-10 386 %% VAL_X00_Y02 F64 8.123932e-08 387 %% VAL_X01_Y02 F64 1.328378e-11 388 %% VAL_X00_Y03 F64 1.170865e-10 389 %% NELEMENTS S32 10 390 %% END 391 %% # These are the original linear solutions 392 %% GHOST.INNER.MAJOR METADATA 393 %% NORDER_X S32 1 394 %% VAL_X00 F64 3.926693e+01 395 %% VAL_X01 F64 5.325759e-03 396 %% NELEMENTS S32 2 397 %% END 398 399 %% GHOST.INNER.MINOR METADATA 400 %% NORDER_X S32 1 401 %% VAL_X00 F64 5.287548e+01 402 %% VAL_X01 F64 -2.191669e-03 403 %% NELEMENTS S32 2 404 %% END 405 406 %% GHOST.OUTER.MAJOR METADATA 407 %% NORDER_X S32 1 408 %% VAL_X00 F64 7.928722e+01 409 %% VAL_X01 F64 1.722181e-02 410 %% NELEMENTS S32 2 411 %% END 412 413 %% GHOST.OUTER.MINOR METADATA 414 %% NORDER_X S32 1 415 %% VAL_X00 F64 1.314265e+02 416 %% VAL_X01 F64 -2.627153e-03 417 %% NELEMENTS S32 2 418 %% END 426 \begin{deluxetable}{lcc} 427 \tablecolumns{3} 428 \tablewidth{0pc} 429 \tablecaption{Optical Ghost Center Transformations} 430 \tablehead{\colhead{Polynomial Term}&\colhead{X center}&\colhead{Y center}} 431 \startdata 432 $x^0 y^0$ & -1.215661e+02 & 2.422174e+01 \\ 433 $x^1 y^0$ & 1.321875e-02 & 4.170486e-04 \\ 434 $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\ 435 $x^3 y^0$ & 1.148288e-10 & -1.173657e-12 \\ 436 $x^0 y^1$ & -1.908074e-03 & 1.189352e-02 \\ 437 $x^1 y^1$ & 8.479150e-08 & -9.256748e-08 \\ 438 $x^2 y^1$ & 1.635732e-11 & 1.140772e-10 \\ 439 $x^0 y^2$ & 2.625405e-08 & 8.123932e-08 \\ 440 $x^1 y^2$ & 1.125586e-10 & 1.328378e-11 \\ 441 $x^0 y^3$ & 2.912432e-12 & 1.170865e-10 \\ 442 \enddata 443 \label{tab:ghost_centers} 444 \end{deluxetable} 445 446 \begin{deluxetable}{lcccc} 447 \tablecolumns{5} 448 \tablewidth{0pc} 449 \tablecaption{Optical Ghost Annulus Axis Length} 450 \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}& \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}} 451 \startdata 452 $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\ 453 $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\ 454 \enddata 455 \label{tab:ghost_radii} 456 \end{deluxetable} 457 458 \begin{figure} 459 \caption{Figure of full FOV showing optical ghosts. Possibly only a few OTAs to illustrate shape deformation.} 460 \end{figure} 419 461 420 462 \subsubsection{Glints} … … 429 471 have a dynamic mask constructed when a reference source falls on the 430 472 focal plane within \czwdraft{approximately} one degree of the detector 431 edge. This mask is 150 pixels wide, and$L = 2500 \left(-20 -432 m_{inst}\right)$. 473 edge. This mask is 150 pixels wide, with length $L = 2500 \left(-20 - 474 m_{inst}\right)$. \czwdraft{Am I correct that this is basically a one-degree edge around the detector?} 433 475 434 476 %% … … 456 498 %% END 457 499 500 \begin{figure} 501 \caption{Example of glint.} 502 \end{figure} 503 458 504 \subsubsection{Diffraction spikes} 459 505 … … 462 508 with length $L = 10^{0.096 * (7.35 - m)} - 200$ and width $W = 8 + (L 463 509 - 200) * 0.01$. These spikes are dependent on the camera rotation, 464 and are oriented at $\theta = n * \frac{ pi}{2} - \mathrm{ROTANGLE} +510 and are oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} + 465 511 0.798$. 466 512 … … 469 515 The cores of saturated stars are masked as well, with radius $r = 10.15 * (-15 - m_{inst})$. \czwdraft{good job here.} 470 516 471 \czwdraft{Write up something about the masking fraction.} 517 \begin{figure} 518 \caption{Example of saturated star, which will also nicely show the diffraction spikes.} 519 \end{figure} 472 520 473 521 \subsection{Video Mask} … … 493 541 \subsection{Masking fraction} 494 542 495 \czwdraft{\% due to chip/cell gaps} 496 497 \czwdraft{\% due to faulty pixels} 498 499 \czwdraft{\% due to CTE} 500 501 \czwdraft{\% due to vinetting} 502 503 \czwdraft{\% average dynamic masking} 543 For the full field of view that falls on the sixty OTAs, 14.7\% \czwdraft{check this} of all pixels are masked. The majority of this masking is due to regions that fall within the vignetted region. Defining the radius of the unvignetted region to be 3 degrees, and excluding pixels that fall beyond this point reduces the static masking fraction to 9.7\%. 544 545 Unfortunately, due to the design of the OTAs and readout cells, a non-negligible fraction of the field of view falls onto an area that does not have a detector pixel. For a given OTA mosaicked to a $4846\times{}4868$ pixel image, the 64 $590\times{}598$ pixel readout cells cover 95.7\% of the OTA area. 546 547 For the inter-chip gap area loss, we use two field of view calculations. The reference field of view of GPC1 is 3 degrees, which at the nominal plate scale of 0.258 arcseconds per pixel, translates to a 20930 FPA pixel radius. However, based on the manual masking of the vignetted region, illuminated pixels are generally unvignetted out to 3.25 degrees, or a 22720 FPA pixel radius. Although these result in different coverage areas, summing the number of pixels on OTA mosaicked images ($4846\times{}4868$ pixels) within either field of view results in a inter-chip gap mask fraction of 7\%. 548 549 %% mysql> select filter,AVG(camProcessedExp.maskfrac_ref_static), AVG(camProcessedExp.maskfrac_ref_dynamic), AVG(camProcessedExp.maskfrac_ref_advisory), AVG(camProcessedExp.maskfrac_max_static),AVG(camProcessedExp.maskfrac_max_dynamic),AVG(camProcessedExp.maskfrac_max_advisory) from camRun join camProcessedExp USING(cam_id) JOIN chipRun USING(chip_id) JOIN rawExp USING(exp_id) WHERE camRun.label = 'LAP.PV3.20140730.final' GROUP BY filter; 550 %% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+ 551 %% | filter | AVG(camProcessedExp.maskfrac_ref_static) | AVG(camProcessedExp.maskfrac_ref_dynamic) | AVG(camProcessedExp.maskfrac_ref_advisory) | AVG(camProcessedExp.maskfrac_max_static) | AVG(camProcessedExp.maskfrac_max_dynamic) | AVG(camProcessedExp.maskfrac_max_advisory) | 552 %% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+ 553 %% static dynamic advisory 554 %% | g.00000 | 0.19642137972007 | 0.00010322263512709 | 0.026838445469766 555 %% | 0.20949461794863 | 9.89200027293e-05 | 0.026431927734548 | 556 %% | r.00000 | 0.19675996201399 | 0.00025214447869606 | 0.032641054600788 557 %% | 0.20989768279138 | 0.00023994155711801 | 0.032178525485201 | 558 %% | i.00000 | 0.19677587604327 | 0.00057470697316504 | 0.038096251937072 559 %% | 0.21003570722292 | 0.00053987093278142 | 0.037471018638997 | 560 %% | z.00000 | 0.1974290315691 | 0.00024758901226967 | 0.03064123748973 561 %% | 0.21055007930696 | 0.00023452690039757 | 0.030144453360769 | 562 %% | y.00000 | 0.19828990634315 | 0.00014523787521897 | 0.021984846417987 563 %% | 0.21130344126869 | 0.00013634812877977 | 0.02163070300815 | 564 565 Summing mask fractions from these three contributions results in an average of $\sim 20\%$ masking fraction across the field of view. Dynamic masking adds an additional $2-3\%$, with advisory burntool masking contributing the largest component. 504 566 505 567 \section{Overscan} … … 524 586 frames with a ramp of exposure times. As the exposure time increases, 525 587 the flux on each pixel also increases in what is expected to be a 526 linear manner. Each of these dark ramp exposuresis overscan527 corrected, and then the median is calculated for each cell, as well as 588 linear manner. Each of these dark exposures in this exposure time ramp is overscan 589 corrected, and then the median is calculated for each cell, as well as for 528 590 the rows and columns within ten pixels of the edge of the science 529 591 region. From these median values at each exposure time value, we can … … 537 599 cell results in the center of the detector. 538 600 539 This non-linearity effect appears to be stable in time, with no540 evident change over a year's worth of data.601 This non-linearity effect appears to be stable in time, with little 602 evident change over the survey duration. 541 603 542 604 \czwdraft{I have figures at http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity that might be useful} 605 %http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity_AllEdges 606 %http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearityArchive 607 608 \begin{figure} 609 \caption{Example plot of linearity as a function of incident brightness.} 610 \end{figure} 543 611 544 612 \section{Dark/Bias Subtraction} 545 613 % http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Background_Dark_Model 546 614 The dark model we make for GPC1 considers each pixel individually, 547 independent of any neighbors. To create the dark model for each 548 pixel, we fit an arbitrary dimensional model \czwdraft{clunky} to the 549 array of input pixels from a selection of dark images. The current 550 model is linear \czwdraft{really?} in both the exposure time and the 551 detector temperature. Adding in a constant value for the fit provides 552 three parameters that define the dark model for that pixel. As this 553 constant value is effectively the bias value for that pixel, we do not 554 do a separate bias correction. This model is applied to science 555 images by fitting the correct dark value based on the exposure time 556 and detector temperature for that exposure. 615 independent of any neighbors. To create the dark model, we fit an multi-dimensional model to the array of input pixels 616 from a randomly selected set of 100-150 \czwdraft{overscan corrected} 617 dark frames chosen from a given date range. The model fits 618 each pixel as a function of the exposure time $t_{exposure}$ and the 619 detector temperature $T_{chip}$ such that $dark = a_0 + a_1 620 t_{exposure} + a_2 T_{chip} t_{exposure} + a_3 T_{chip}^2 621 t_{exposure}$. This fitting is performed over the sample of input pixels, 622 and the coefficients $a_i$ stored in the detrend image. The constant 623 $a_0$ term includes the bias signal, and as such, a separate bias 624 subtraction is not necessary. 625 626 Applying the dark model is simply a matter of calculating the response 627 for the exposure time and detector temperature for the image to be 628 corrected, and subtracting the resulting dark signal from the image. 557 629 558 630 \subsection{Time evolution} … … 560 632 \czwdraft{The dark model is noticably unstable on time scales of months, and so we have generated a sequence in time to keep the effect of a missed correction low.} 561 633 562 Unfortunately, the dark model is not consistently stable on the time 563 span of multiple months. Some of the changes in the dark can be 634 The dark model is not consistently stable over the full survey, with significant drift over the course of multiple months. Some of the changes in the dark can be 564 635 attributed to changes in the voltage settings of GPC1, but the 565 majority seem to be the result of some unknown parameter. Largely, we636 majority seem to be the result of some unknown parameter. We 566 637 can separate the dark model history of GPC1 into three epochs. The 567 638 first epoch covers all data taken prior to 2010-01-23. This epoch … … 572 643 characterized by a largely stable but oscillatory dark solution. 573 644 There appear to be two modes that the dark model switches between 574 apparently at random. No clear cause has been established for the se645 apparently at random. No clear cause has been established for the 575 646 switching, but there are clear differences between the two modes 576 647 \czwdraft{figures?}. … … 586 657 appropriate ``A'' and ``B'' mode dark frames. Using the appropriate 587 658 dark minimizes the effect of this bias gradient in the dark corrected 588 data. 659 data. Table \ref{tab:dark mode dates} lists the dates used for each dark mode. 660 661 \czwdraft{The evidence of the mode switching can be visualized in Figure \ref{fig:dark switching}. This figure shows image profile along the x-pixel axis binned along the full y-axis of dark corrected images for OTA67. These images are from sequential days, and have been corrected with a dark model constructed from the full set of dark data within the second epoch. The opposite sign of the slopes of these profiles indicates that the average dark model does not correct these dates sufficiently, due to the contradictory dark signals between the two modes.} 589 662 590 663 After 2011-05-01, the two-mode behavior of the dark disappears, and is … … 594 667 These darks cover the range from 2011-05-01 to 2011-08-01, 2011-08-01 595 668 to 2011-11-01, and 2011-11-01 and on. The reason for this time 596 dependent driftis unknown, but we seem to be able to model it with669 evolution is unknown, but we seem to be able to model it with 597 670 reasonable accuracy by creating new dark models. 598 671 672 \begin{figure} 673 \caption{Example of raw and dark calibrated exposure. Plots of horizontal cuts for A/B/average corrections.} 674 \end{figure} 675 676 \subsection{Video Dark} 677 678 Dark signal is stronger in cell corners due to amplifier glow. Standard model corrects this. When OTA reads video cell, the dark model changes. The standard model oversubtracts the dark model. Make video darks from dark data that has had video signal running. Need two passes to cover all cells (shifting video cell between the two). Can construct the video dark and the standard dark simultaneously, by using OTAs that have video on and off. 679 680 Video dark signal appears linear and stable, allowing archival data from prior to video dark data to be corrected by simply taking $VD_{2009} = D_{2009} - D_{Modern} + VD{Modern}$. 681 599 682 \section{Noisemap} 600 683 601 Based on a study of the positional dependence of detected objects, we discovered that the cells in GPC1 do not have uniform noise characteristics. Instead, there is a gradient along the pixel rows, with the noise generally higher away from the read out amplifier. This is likely another effect of the row-by-row bias issue. This gradient has the effect that the read noise increases as the row is readout. To mitigate this noise gradient, we construct a set of noisemap images by measuring the median variance on bias frames. The variance is calculated in boxes of 20x20 pixels, and then linearly interpolated to cover the full image.602 603 Unfortunately, due to correlations in the row-to-row offsets \czwdraft{in the noise?}, the variance measured from the bias images does not fully remove the positional dependence of objects that are detected. The reason for this is that the simple noisemap underestimates t ehnoise observed when the image is filtered during the object detection process. This filtering convolves the background noise with a PSF, which has the effect of amplifying the correlated peaks in the noise. This amplification can therefore boost background fluctuations above the threshold used to select real objects, contaminating the final object catalogs.604 605 To resolve this issue, we chose a typical PSF, and used it to look for detections on a sample of bias images. As the bias has no real sources, all objects found are by definition false, and provides an idea of how much our noisemap estimation deviates from the ``true'' noise observed by the object detection process. For a region of area X*Y, if we find k false detections above our signal-to-noise threshold, then we can estimate how much the noise model deviates from what is observed. The observed noise threshold is defined as $\sigma_{observed} = \sqrt{2} * \erfcinv{2 * k A_{psf} / (X * Y * N_{exp})}$, where $A_{psf}$ is the footprint size of the PSF (taken as 16 pixels), and $N_{exp}$ is the number of exposures examined in this location. From this observed threshold, we scale the noisemap previously calculated by the boost factor $B = \sigma_{thresh} / \sigma_{observed}$.606 607 The row-to-row variations that contribute to the extra noise are related to the dark model, and because of this, as the dark model changes, the effective noise also changes. Because of this, we have created different noisemap models for the three major time ranges of the dark model. We do not see any evidence that the noisemaps have the A/B modes visible in the dark, and so we do not generate different models.684 Based on a study of the positional dependence of detected sources, we have discovered that the cells in GPC1 do not have uniform noise characteristics. Instead, there is a gradient along the pixel rows, with the noise generally higher away from the read out amplifier. This is likely another effect of the row-by-row bias issue discussed below. This gradient has the effect that the read noise increases as the row is read out. As a result of this increased noise, more sources are detected when the readnoise is assumed constant across the readout. To mitigate this noise gradient, we construct a set of noisemap images by measuring the median variance on bias frames. The variance is calculated in boxes of 20x20 pixels, and then linearly interpolated to cover the full image. 685 686 Unfortunately, due to correlations in the row-to-row offsets \czwdraft{in the noise?}, the variance measured from the bias images does not fully remove the positional dependence of objects that are detected. The reason for this is that the simple noisemap underestimates the noise observed when the image is filtered during the object detection process. This filtering convolves the background noise with a PSF, which has the effect of amplifying the correlated peaks in the noise. This amplification can therefore boost background fluctuations above the threshold used to select real objects, contaminating the final object catalogs. 687 688 To resolve this issue, we chose a PSF with a typical FWHM, and used it to look for detections on a sample of bias images. As the bias has no real sources, all objects found are by definition false, which provides an idea of how much our noisemap estimation deviates from the ``true'' noise observed by the object detection process. For a region of area $X*Y$, if we find $k$ false detections above our signal-to-noise threshold $sigma_{thresh}$, then we can estimate how much the noise model deviates from what is observed. The observed noise threshold is defined as $\sigma_{observed} = \sqrt{2} * \erfcinv{2 * k A_{psf} / (X * Y * N_{exp})}$, where $A_{psf}$ is the footprint size of the PSF (taken as 16 pixels), and $N_{exp}$ is the number of exposures examined in this location. From this observed threshold, we scale the noisemap previously calculated by the boost factor $B = \sigma_{thresh} / \sigma_{observed}$. 689 690 The row-to-row variations that contribute to the extra noise are related to the dark model, and because of this, as the dark model changes, the effective noise also changes. To ensure that the noisemap accurately matches the true noise level, we have created different noisemap models for the three major time ranges of the dark model. We do not see any evidence that the noisemaps have the A/B modes visible in the dark, and so we do not generate different models for each individual dark model. 608 691 609 692 \section{Remnance?} … … 623 706 uniformly illuminated image. Using a dome screen is not possible, as 624 707 the variations in illumination and screen rigidity create unusably 625 large scatter between different images . Because of this, we use sky708 large scatter between different images that are caused by the detector response function. Because of this, we use sky 626 709 flat images taken at twilight, which are more consistently illuminated 627 710 than screen flats. We calculate the mean of these images to determine … … 630 713 From this initial flat model, we construct a correction to remove the 631 714 effect of the problems illuminating the large area. This is done by 632 dithering a series of exposures across a given pointing. By comparing 633 the measured fluxes for a given star as a function of position, we can 634 correct out the errors in the flat model. 715 dithering a series of science exposures across a given pointing. By 716 comparing the measured fluxes for a given star as a function of 717 position on the detector, we can determine the position dependent 718 scaling factors. These scale factors can then be used to correct the 719 initial flat field model to better represent the detector response. 635 720 636 721 The flat model appears stable with time, although directly measuring 637 722 this is as difficult as originally constructing the model. However, 638 due to the photometric consistency observed in GPC1 measurements, we 639 can be confident that the flat model is not changing much. 723 due to the photometric consistency observed in the catalog of GPC1 measurements, we 724 can be confident that the flat model is not as time dependent as the 725 dark correction. 640 726 641 727 … … 647 733 to remove the remaining row-by-row variation, and the PATTERN.CELL and 648 734 PATTERN.CONTINUITY corrections attempt to ensure that the cells of a 649 given OTA are consistent with each other. These corrections are735 given OTA are consistent with the other cells on that OTA. These corrections are 650 736 largely designed to fix issues that are not stable enough with time 651 737 for the dark model or flat field model to fully account for the … … 653 739 654 740 \subsection{Pattern Row} 655 741 % http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study 656 742 As discussed above in the dark and noisemap sections, certain 657 detectors have significant row-by-row bias offsets. As the level of 658 the offset is largely random, the dark correction cannot fully remove 659 this structure from the images. Therefore, we apply the PATTERN.ROW 660 correction in an attempt to mitigate the offsets. To force the rows 661 to agree, a \czwdraft{first} order polynomial is fit to each row in the 662 cell, and that trend subtracted from the data. The median offset 663 (corresponding to the background level) is then added back to the 664 image so that the cell matches its neighbors during background 665 subtraction. 743 detectors have significant row-by-row bias offsets. The magnitude of 744 these offsets increases as the distance from the readout amplifier 745 increases, resulting in horizontal streaks that are more pronounced 746 along one edge of the cell. As the level of the offset is largely 747 random, the dark correction cannot fully remove this structure from 748 the images, and the noisemap level only indicates the level of the 749 variance added by these bias offsets. Therefore, we apply the 750 PATTERN.ROW correction in an attempt to mitigate the offsets. To 751 force the rows to agree, a \czwdraft{second} order polynomial is fit to 752 each row in the cell, and that trend subtracted from the data. The 753 median offset (corresponding to the background level) for each row is then fit by a first order polynomial, and that trend is then added 754 back to the image so that the sky level on the cell matches its neighbors during 755 background subtraction. 666 756 667 757 This correction was required on all cells on all OTAs prior to … … 669 759 electronics resolved the row-by-row offsets for the majority of the 670 760 detectors. As a result, we only apply this correction where it is 671 necessary, as shown in figure \czwdraft{X}.761 necessary, as shown in Figure \ref{fig: pattern row required}. 672 762 673 763 Although this correction does resolve the row-by-row offset issue in a … … 676 766 near these objects. As the offsets are calculated on the pixel rows, 677 767 this oversubtraction is not uniform around the object, but is 678 preferentially along the $\pm x$ axis of the object. 679 680 \czwdraft{keep this?} This row-by-row offset is visible in similar 681 camera designs, and has been removed by identifying the noise signal 682 in the pixel data stream. By taking the FFT of the pixels and a 683 reference signal, the frequency of this noise can be isolated and 684 removed, resulting in a much cleaner image. However, GPC1 does not 685 record the value of the reference signal, instead automatically 686 subtracting it from the data values. Without this comparison signal, 687 we have been unable to reproduce this method, as there is no obvious 688 FFT component visible. 768 preferentially along the horizontal x axis of the object. 769 770 %% \czwdraft{keep this?} This row-by-row offset is visible in similar 771 %% camera designs, and has been removed by identifying the noise signal 772 %% in the pixel data stream. By taking the FFT of the pixels and a 773 %% reference signal, the frequency of this noise can be isolated and 774 %% removed, resulting in a much cleaner image. However, GPC1 does not 775 %% record the value of the reference signal, instead automatically 776 %% subtracting it from the data values. Without this comparison signal, 777 %% we have been unable to reproduce this method, as there is no obvious 778 %% FFT component visible. 779 780 \begin{figure} 781 \caption{Example of pre/post pattern row application.} 782 \end{figure} 689 783 690 784 \subsection{Pattern Cell} … … 692 786 As the bias level of a given cell may not exactly match that of its 693 787 neighbors, fitting a smooth background model results in over and 694 under-subtraction of the sky level at these discontinuities. The 695 PATTERN.CELL correction was the first attempt to remove this effect on 696 the worst cells, by forcing all the cells of an OTA to the same level. 697 Each cell has the median value measured, and then an offset added that 698 shifts each cell to match the median of these medians. 788 under-subtraction of the sky level at the cell boundary 789 discontinuities. The PATTERN.CELL correction was the first attempt to 790 remove this effect on the worst cells, by forcing all the cells of an 791 OTA to the same level. Each cell has the median value measured, and 792 then each cell has an offset added that shifts the cell to match the 793 median of those medians. 699 794 700 795 This correction is reasonable when the astronomical signal is smooth, 701 796 with no objects that are large relative to the size of an individual 702 797 cell. However, the presence of large galaxies (or even bright stars) 703 can force some cells into a nearly arbitrary offset from their 704 neighbors. Because of this issue, we no longer apply this correction 705 to any data. 798 can bias the offsets for some cells from their neighbors. Because of 799 this issue, we no longer apply this correction to any data. 706 800 707 801 \subsection{Pattern Continuity} 708 802 709 As the PATTERN.CELL correction was clearly defectivein many710 situations, we designed a replacement correction that would distort711 large objects less. In addition, studies of the background level 712 illustrated that the row-by-row bias introduces a small background 713 gradient along the rows of the cells. This results in a ``sawtooth'' 714 pattern across an OTA, and as the background model assumes a smooth 715 sky level, we saw evidence of over and under subtraction at cell 716 boundaries. As the PATTERN.CELL was designed to correct mean changes 717 between cells, it could not adequately resolve this higher order 718 issue.803 As the PATTERN.CELL correction was clearly insufficient in many 804 situations, we designed a replacement correction that would lower the 805 distortion for large objects less. In addition, studies of the 806 background level illustrated that the row-by-row bias introduces 807 small background gradient variations along the rows of the cells that is not stable enough to be completely fit by the dark model. This results 808 in a ``sawtooth'' pattern horizontally across an OTA, and as the background model 809 assumes a smooth sky level, this induces over and under 810 subtraction at cell boundaries. As the PATTERN.CELL was designed to 811 correct mean changes between cells, it could not adequately resolve 812 this higher order issue. 719 813 720 814 The replacment for PATTERN.CELL was the PATTERN.CONTINUITY correction, 721 815 which attempts to match the edges of a cell to those of its neighbors. 722 For each cell, a thin box on each edge is extracted and the median816 For each cell, a thin box \czwdraft{10} pixels wide on each edge is extracted and the median 723 817 value calculated for that box. These median values are then used to 724 818 construct a vector of differences $diff_i = \sum_{j,j'} Edge_{i,j} - 725 819 Edge)_{i',j'}$, along with a matrix of associations $A_{i,i'} = 726 \sum_{j,j'} \delta(j,j')$ denoting which cell boundar y touches820 \sum_{j,j'} \delta(j,j')$ denoting which cell boundaries touch 727 821 another. By solving the system $A x = diff$, we can find the set of 728 822 offsets $x_i$ that should be applied to each cell to ensure the … … 732 826 to align the cells into a single ramp, at the expense of the absolute 733 827 background level. However, as we subtract off a smooth background 734 model, th is absolute level isunimportant. The fact that the final828 model, the deviations from an absolute sky level are unimportant. The fact that the final 735 829 ramp is smoother than it would be otherwise also allows for the 736 830 background subtracted image to more closely match the astronomical 737 831 sky, without over- and under-subtractions at cell edges. 738 832 739 %% \section{Fringe correction} 740 741 %% \czwdraft{Due to variations in the thickness of the detectors, we observe interference patterns at the infrared (red?) end of the filters, as the wavelength of the light becomes comparable to these variations. Visually inspecting the images shows that the fringing is most prevalent in the y-filter images, with minimal fringing in other bands. Because of this, we only apply a fringe correction to the y data.} 742 743 %% \czwdraft{The fringe is constructed by randomly determining a set of boxes for each OTA cell, and measuring the sky subtracted median value in those boxes for a series of images. These samples are selected at the same location on each image, allowing the astronomical signal to be removed. A least squares fit to the data is then calculated, providing the model of the fringe strength at that location.} 744 745 %% \czwdraft{Applying the fringe is done in the same way, with samples measured across the image to determine the relative strength of the fringing in this image. The solution derived from the detrend is then scaled to match that observed in the science image, and subtracted away.} 746 747 %% \section{Background subtraction} 748 749 %% \czwdraft{A background model is generated for each OTA, once all the individual cells have been mosaicked together. Super-pixels are then defined that divide the image into XxY subregions, and the mean calculated for each subregion. This grid is shifted by a half-width, and the means recalculated, to double the sampling frequency. A background model is then calculated by interpolating over this sampled grid.} 833 \begin{figure} 834 \caption{Continuity example, with background issue.} 835 \end{figure} 836 837 \section{Fringe correction} 838 839 Due to variations in the thickness of the detectors, we observe 840 interference patterns at the infrared end of the filter set, as 841 the wavelength of the light becomes comparable to the thickness of these variations. 842 Visually inspecting the images shows that the fringing is most 843 prevalent in the y-filter images, with minimal fringing in other 844 bands. As a result of this, we only apply a fringe correction to the y filter 845 data. 846 847 The fringe is constructed by randomly determining a set of boxes for 848 each OTA cell, and measuring the sky subtracted median value in those 849 boxes for a series of images. These samples are selected at the same 850 location on each image, allowing the astronomical signal to be 851 filtered out as an additional noise term. A least squares fit to the 852 data is then calculated, providing the model of the fringe strength at 853 that location. 854 855 Applying the fringe is done in the same way, with samples measured 856 across the image to determine the relative strength of the fringing in 857 this image. The solution derived from the detrend is then scaled to 858 match that observed in the science image, and subtracted away. 859 860 \section{Background subtraction} 861 862 \czwdraft{A background model is generated for each OTA, once all the individual cells have been mosaicked together. Super-pixels are then defined that divide the image into XxY subregions, and the mean calculated for each subregion. This grid is shifted by a half-width, and the means recalculated, to double the sampling frequency. A background model is then calculated by interpolating over this sampled grid.} 863 864 %% * Magic 865 %% * Warping 866 %% * warping kernel 867 %% * linear-by-pieces 868 %% * Covariance 869 %% * def of skycells? 870 %% * Stacking 871 %% * pixel combination rules 872 %% * pixel rejections 873 %% * convolution for matching (success and failure) 874 %% * Difference Image analysis 875 876 877 \section{Warping} 878 879 To provide a consistent and uniform set of images for co-added image 880 stacking and image differences, the individual mosaicked OTA images 881 are projected onto a common set of tangent plane projected regions. 882 These projection cells are $4\times{}4$ degree fields spaced onto set 883 of projection centers that fully cover the sky. These projection 884 cells are arranged into rings of constant declination, and allowed to 885 overlap as $|\delta|$ increases. Each projection cell is further 886 subdivided into \czwdraft{size} sky cells, which have constant overlap 887 regions of \czwdraft{overlap}. These skycells are the main image unit 888 used for processing image data beyond the initial chip stage. 889 890 After the detrending and photometry, the detection catalog for the 891 full camera is fit to the reference catalog, producing third-order 892 astrometric solutions that map the detector focal plane to the sky, 893 and map the individual OTA pixels to the detector focal plane. This 894 solution is then used to determine which skycells the exposure OTAs 895 overlap. 896 897 Foreach output skycell, all overlapping OTAs and the calibrated 898 catalog are read into the \textbf{pswarp} program. Each input image 899 is examined in order, and the same transformation performed. This 900 transformation breaks the output warp image into $128\times{}128$ 901 pixel grid boxes. Each grid box has a locally linear map calculated 902 that converts the output warp image coordinates to the input chip 903 image coordinates. By doing the transformation in this direction, 904 each output pixel has a unique sampling position on the input image 905 (although it may be off the image frame and therefore not populated), 906 preventing gaps in the output image due to the spacing of the input 907 pixels. 908 909 With the locally linear grid defined, Lanczos interpolation with 910 filter size parameter $a = 3$ on the input image is used to determine 911 the values to assign to the output pixel location. The output 912 locations are shifted by 0.5 pixels to let the interpolation select 913 the value that would be assigned to the center of the output 914 pixel. This process is repeated for all grid boxes, for all input 915 images, and for each output image product: the science image, the 916 variance, and the mask. \czwdraft{The jacobian is multiplied to the 917 image value, and squared and multiplied to the variance. I don't 918 understand that.} 919 920 As the interpolation constructs the output pixels from more than one 921 input pixel, there is a covariance term that is must be included. For 922 each locally linear grid box, the covariance is calculated from the 923 kernel in the center of the 128 pixel range. Once the image has been 924 fully populated, this set of individual covariance matrices is 925 averaged to create the final covariance for the full image. 926 927 An output catalog is also constructed from the full exposure input 928 catalog, including only those objects that fall on the warped image. 929 These detections are transformed to match the new image location, and 930 to scale the position errors based on the new orientation. 931 932 The output image also contains header keywords SRC\_0000, SEC\_0000, 933 MPX\_0000, and MPY\_0000 that contain the mappings from the warped 934 pixel space to the input image. The SRC keyword lists the input OTA 935 name, and the SEC keyword lists the image section corresponding to the 936 locally linear grid box. The MPX and MPY contain the transformation 937 parameters for the locally linear grid. \czwdraft{Is this accurate?} 938 939 % Read all images and astrometry 940 % Check which input images overlap with output image. => 8007 when the inputs don't overlap. 941 % Loop over each image. 942 % Append detections from input into output detection list. 943 % Determine transform back from warp pixels to source image. 944 %% 2nd order polynomial in both x and y for this transformation. and save to header 945 % Break warp image into 128x128pixel locally linear areas 946 % Mask non finite pixels as saturated. 947 % Define Lanczos-3 interpolation over the input image. 948 % Iterate over warp pixel space (on each locally linear grid) and map interpolated input pixel positions onto warp. 949 % Warp pixel space is defined as center based, so that's where the intpolation point comes from. 950 % Covariance calculated based on the interpolation kernel at the center of the ll grid. 951 % image = interp_image * jacobian 952 % var = interp_var * jacobian**2 953 % mask = interp_mask 954 % jacobian = abs(mapXx * mapYy - mapYx * mapXy) 955 % I don't understand that jacobian. 956 % 957 958 959 \section{Stacking} 960 961 Once individual exposures have been warped onto a common projection 962 system, they can then be combined without that added concern. In 963 order to obtain detections of faint images, and to provide a static 964 sky image without transient features, we coadd the individual warps 965 into a stacked image. Creating this stack also allows a complete 966 image to be constructed that does not have regions masked due to 967 falling between devices. 968 969 The stacked image is comprised of all warp frames for a given skycell 970 in a single filter. The source catalogs and images are loaded into 971 the \textbf{ppStack} program to do prepare the inputs and stack the 972 frames while rejecting bad pixels. 973 974 Once all files are ingested, the first step is to measure the size and 975 shapes of the input image PSFs. We exclude images that have a PSF 976 FWHM greater than 10 pixels, as those images have the worst seeing and 977 would degrade the final output stack. A target PSF for the stack is 978 constructed from the envelope of all input PSFs, which sets the target 979 PSF at the largest value among the input PSFs for all radii. This PSF 980 is then circularized to prevent any of the input images from being 981 deconvolved when matched to the target. 982 983 The input images also need to be normalized to prevent differences in 984 seeing and sky transparency from causing discrepancies during pixel 985 rejection. From the calibrated input catalogs, we have the 986 instrumental magnitudes of all sources, along with the airmass, image 987 exposure time, and zeropoint. All output stacks are calibrated to a 988 zeropoint of 25.0 in all filters, and to have an airmass of 1.0. The 989 output exposure time is set as the sum of the input exposure times. 990 With this information, we can determine the relative transparency for 991 each input image by comparing matched sources between the different 992 images. Each image then has a normalization factor defined, equal to 993 $norm_{image} = (ZP_{image} - ZP_{target}) - transparency_{image} - 994 2.5 * \log_{10} (t_{target} / t_{image}) - airmassTerm * 995 (airmass_{image} - airmass_{target})$. The input source catalog is 996 adjusted to reflect this normalization, which is also retained for 997 application when the pixels are combined. 998 999 % PREPARE 1000 % load sources 1001 % load psfs 1002 % determine target as envelope of input psfs 1003 % FWHM clipping at 10 1004 % measure seeing 1005 % - // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency 1006 % // EAM : the discussion here was not quite right (or at least sloppy). Here is a replacement explanation: 1007 % // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by: 1008 % // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency 1009 % // NOTE the sign of 'transparency' this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where 1010 % // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app 1011 % // we want to adjust the input images to be in a consistent flux system so that the 1012 % // final stack can be generated with a specific target zero point. Any adjustment to 1013 % // the flux scale of the image must be made in coordination with the resulting 1014 % // zeropoint, exposure time, and airmass such that the above relationship yields the 1015 % // same apparent magnitude for a given star: 1016 % // m_inst_i : instrumental mags on input image (in) 1017 % // m_inst_o : instrumental mags on re-normalized image (out) 1018 % // m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i 1019 % // m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o 1020 % // zp_i, airmass_i, t_i, trans_i : reported or measured for input image 1021 % // zp_o = zpTarget (from recipe) 1022 % // airmass_o = airmassTarget (from recipe) 1023 % // t_o = sumExpTime [sum of input exposure times: once images are scale to this time, they can be avereaged] 1024 % // trans_o = 0.0 [obviously!] 1025 % // we have 2 cases: (a) all reported ZPs are good or (b) some are bad: 1026 % // (a) FPA.ZP = zp_i + c1 * airmass_i 1027 % // --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i) 1028 % // (b) zp[i] = c1 * airmass_i + 2.5log(exptime_i) 1029 % // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong. 1030 % // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i 1031 1032 With the normalization factors and target PSF chosen, the convolution 1033 kernels can be calculated for each image. ISIS kernels are used with 1034 FWHM values of 1.5, 3.0, and 6.0 pixels and polynomial orders of 6, 4, 1035 and 2. \czwdraft{Skipping this bit because I'm not completely sure I 1036 understand it.} The image is then scaled by the normalization as 1037 $renorm = 10^{-0.4 * norm_{image}} / norm_{convolution}$, and the 1038 variance by the square of that value. 1039 1040 % MATCH 1041 % match to target PSF. 1042 % use ISIS kernels to do matching/convolution 1043 % Input sources used for psf matching. 1044 % @ISIS.WIDTHS F32 1.5 3.0 6.0 # Gaussian kernel FWHM values 1045 % @ISIS.ORDERS S32 6 4 2 # Polynomial orders for ISIS kernels 1046 1047 Once the convolution kernels are defind for each image, they are used 1048 to convolve the image to match the target PSF. Any input image that 1049 has a $\chi^2$ value larger than 4.0$\sigma$ larger than the median 1050 value is rejected from the stack. 1051 1052 % CONVOLVE 1053 % Normalization to match target zeropoint/exptime 1054 % Reject images with bad match chi^2 values. MATCH.REJ * rms + median threshold. 1055 % Additional variance from the convolution chi^2 1056 % Calculate image weights based on variance: W_i = 1 / (ROBUST_MEDIAN(variance_image_i) * CovarianceFactor) 1057 % CovarianceFactor = covariance->kernel[0][0] 1058 1059 Following the convolution, and initial stack is constructed. For a 1060 given pixel coordinate, the values at that coordinate are extracted 1061 from all input images. Images that have a suspect mask bit (including 1062 the SUSPECT, BURNTOOL, SPIKE, STREAK, STARCORE, and CONV.POOR bit 1063 values) are appended to a suspect pixel list for preferential 1064 exclusion. Following this, the pixel values are combined and tested 1065 to attempt to identify discrepant values that should be excluded. 1066 1067 If only a single input is available, the initial stack contains the 1068 value from that single input. If there are only two inputs, the 1069 average of the two is used. These cases should occur only rarely in 1070 the $3\Pi$ survey, as there are many input exposures that overlap any 1071 particular point on the sky. The more common case for three or more 1072 inputs constructs a weighted average from the inputs, with the weight 1073 set as a single value for each input image, and defined as the inverse 1074 of the median variance value from that image's associated variance 1075 map. This weight is used for the image and the exposure weighted 1076 image: 1077 1078 \begin{eqnarray} 1079 S_{value} &=& \sum_i\left(value_{i} * weight_i\right) / \sum_i\left(weight_i\right) \\ 1080 S_{exp weight} &=& \sum_i \left(exptime_i * weight_i\right) / \sum_i\left(weight_i\right) \\ 1081 \end{eqnarray} 1082 1083 The pixel exposure time is simply the sum of the input exposure time values, and the output variance is 1084 1085 \begin{eqnarray} 1086 S_{variance} &=& 1 / \sum_i \left( 1 / variance_i \right) 1087 \end{eqnarray} 1088 1089 The output mask value is taken to be zero (no masked bits), unless 1090 there were no valid inputs, in which case the BLANK mask bit is set. 1091 1092 % INITIAL COMBINE 1093 % Calculate weighted mean of input images 1094 % mu = sum_i(f_i * W_i) / sum_i(W_i) 1095 % sigma = 1 / sum_i(1 / W_i) 1096 % nu = sum_i(m_i) 1097 % // We're not using the input pixel variances to generate a weighted average for the pixel flux (because 1098 % // that introduces systematic biases), so the variance of the output pixel value should simply be: 1099 % // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 1100 % // This reduces, when the weights are all identically unity, to: 1101 % // variance_combination = sum(variance_i) / N^2 1102 % // and if the variances are all equal: 1103 % // variance_combination = variance_individual / N 1104 % // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). 1105 % sumValueWeight = sum_i(values * weights) 1106 % sumWeight = sum_i(weights) 1107 % sumVarianceWeight == sum( 1 / variances) 1108 % sumExp = sum_i(exptimes) 1109 % sumExpWeight = sum_i(exptims * weights) 1110 % mean = sumValueWeight / sumWeight 1111 % var = 1 / sumVarianceWeight 1112 % exp = sumExp 1113 % expWeight = sumExpWeight 1114 1115 % EXCEPT: if N = 1, accept it. if N = 2, take average. 1116 1117 % if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter, 1118 % combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) { 1119 %bool pmStackCombine( 1120 % pmReadout *combined, // output stacked readout 1121 % pmReadout *expmaps, // output exposure map information 1122 % psArray *input, // input exposures 1123 % psImageMaskType badMaskBits, // treat these bits as 'bad' 1124 % psImageMaskType suspectMaskBits, // treat these bits as 'suspect' 1125 % psImageMaskType blankMaskBits, // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) 1126 % int kernelSize, 1127 % float iter, 0.5 1128 % float rej, 4.0 1129 % float sys, 0.1 1130 % float olympic, 0.2 1131 % bool useVariance, 1132 % bool safe, 1133 % int nminpix, 1134 % bool rejection) 1135 %{ 1136 1137 % combineExtract 1138 %% pixels with mask values as suspect are appended to suspect pixel list. 1139 % combinePixels 1140 %% As described above. 1141 1142 Following this initial combination, a ``testing'' loop iterates in an 1143 attempt to identify outlier points. Again, if only one input is 1144 available, that input is accepted. If there are two inputs, $A$ and 1145 $B$, then a check is made to see if $(0.5 * (value_A - value_B))**2 > 1146 rej**2 * (variance_A + variance_B + (sys * value_A)**2 + (sys * 1147 value_B)**2)$, where $rej$ is the number of sigma deviant a point 1148 needs to be to be excluded, set to 4.0 for the PV3 processing, and 1149 $sys$ is an estimate of the systematic error, taken to be 0.1. 1150 1151 1152 \czwdraft{This discussion seems out of place, but I'm not sure where a 1153 better place is.} Due to the various non-astronomical ghosts that 1154 can occur on GPC1, and the fact that they may not be masked 1155 aggressively enough to ensure all bad pixels are removed, it is 1156 expected that some of the inputs for a given stack pixel are not in 1157 agreement with the others. In general, there is the population of 1158 input pixel values around the correct astronomical level, as well as 1159 possible populations at lower pixel value (such as due to an 1160 over-subtracted burntool trail) and at higher pixel values (such as 1161 that caused by an incompletely masked optical ghost). Due to the 1162 observation strategy to image a given field twice to allow for 1163 warp-warp difference images to be constructed to identify transient 1164 detections, higher pixel values that come from sources like optical 1165 ghosts that are a function of pointing will come in pairs as well. 1166 The higher pixel value contaminants are also potentially problematic 1167 as they may appear to be a real source, prompting photometry to be 1168 performed on a false object. Because of these reasons, there is a 1169 slight preference to reject higher pixel values than lower pixel 1170 values. 1171 1172 If the number of inputs is larger than 6, then a Gaussian mixture 1173 model analysis is run on the inputs to fit two sub populations, and 1174 determine an the likelihood that the distribution is best described by 1175 an uni-modal model. If this probability is less than 0.05, then the 1176 mean is taken from the bimodal sub population with the largest 1177 fraction of inputs, as this should exclude any sub population 1178 comprised of high pixel value outliers. 1179 1180 If this is not the case (the distribution is likely unimodal) or if 1181 there are insufficient inputs for the mixture model analysis, the 1182 input values are passed to an Olympic weighted mean calculation. We 1183 set 0.2 as the fraction of the number of inputs to reject through this 1184 process. This sets the number of bad inputs at $N_{bad} = 0.2 * 1185 N_{input} + 0.5$, where the 0.5 term ensures at least one input is 1186 rejected. This number is further separated into the number of low 1187 values to exclude $N_{low} = N_{bad} / 2$, which will default to zero 1188 if there are few inputs due to integer arithmatic, and $N_{high} = 1189 N_{input} + N_{low} - N_{bad}$. After sorting the input values to 1190 determine which values fall into the low and high groups, the 1191 remaining input values have a weighted mean calculated as described 1192 above. 1193 1194 A systematic variance term is necessary to correctly scale how 1195 discrepant points can be from the ensemble mean. If the mixture model 1196 analysis was run, the Gaussian sigma from the largest sub population 1197 is squared and used. If this is not available, a 0.1 scaling on the 1198 input values is used. Each point then has a limit calculated: 1199 1200 \begin{eqnarray} 1201 limit_{mixture_model} &=& rej**2 * (variance_i + \sigma_{MM}^2) \\ 1202 limit_{default} &=& rej**2 * (variance_i + (0.1 * value_i)**2) 1203 \end{eqnarray} 1204 1205 where $rej$ is the same factor of 4.0 used above. Each input pixel is 1206 then compared against this limit, and the most discrepant pixel that 1207 has $(value_i - mean)**2$ exceeding this limit is identified. If 1208 there are suspect pixels in the set those pixels are marked for 1209 rejection, otherwise this worst pixel is marked for rejection. 1210 Following this, the combine and test loop is repeated for a total $0.5 1211 N_{input}$ iterations, or until no more pixels are rejected. 1212 1213 % combineTest 1214 %% if (Ninput > 6) { use KMM } 1215 %% KMM: 1216 %% Calculate KMMmu KMMsigma KMMpi KMMPunimodal 1217 %% SumWeights = sum(pixelWeights) 1218 %% SysVar = KMMSigma**2 OR (sys * pixelData[i])**2 1219 %% pixelLimts[i] = rej**2 * (pixelVariances[i] + sysVar) 1220 % Iterate 0.5 * Ninput times (at least once) 1221 %% Ninput = 1 => accept 1222 %% Ninput = 2 => if (0.5 * (A - B))**2 > rej**2 * (pixelVariance[A] + pixelVariance[B] + (sys * A)**2 + (sys * B)**2) 1223 %% then if (suspect) mark reject else mark inspect 1224 %% Else => if (useKMM and Punimodal < 0.05) median = KMMmean 1225 %% => else median = combinationWeightedOlympic{} 1226 %% => if (pixelData - median)**2 > pixelLimits[i] then find single worst deviant pixel value 1227 %% then => if suspect (mark reject) else (mark reject worst deviant pixel value) 1228 1229 1230 %% combinationWeightedOlympic => 1231 %% numBad = frac * Ninput + 0.5 1232 %% low = numBad / 2, high = low + numGood - numBad 1233 %% sort(values) => 1234 %% if (i > low && i <= high) { sumValues = sum_i(values * weights); sumWeight = sum_i(weights) 1235 %% return (sumValues / sumWeight) 1236 1237 % obtain lists of inspect and reject pixels. 1238 1239 % normalize:? 1240 % float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation 1241 % psBinaryOp(ro->image, ro->image, ``*'', psScalarAlloc(normalise, PS_TYPE_F32)); 1242 % psBinaryOp(ro->variance, ro->variance, ``*'', psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32)); 1243 1244 With the initial list of rejected pixels generated, a rejection mask 1245 is made by constructing an empty image that has the rejected pixels 1246 set to a value of 1.0. This image is then convolved with a 5 pixel 1247 FWHM 0-order ISIS kernel. Any pixels that are above the threshold of 1248 0.5 are marked as bad and will be rejected in the final convolution. 1249 If more than 10\% of all pixels from an input image are rejected, then 1250 that entire image is rejected as well. 1251 1252 % PIXEL REJECTION 1253 % Construct 15-pixel wide ISIS kernel with 5 pixel FWHM 0-order. 1254 % Construct image of pixels to inspect and convolve with kernel (normalize out kernel power) 1255 % Determine pixels are bad if they're larger than THRESHOLD.MASK = 0.5. 1256 % If more than IMAGE.REJ = 0.1 fraction of pixels are rejected, the entire image is rejected. 1257 1258 1259 \czwdraft{I'm not entirely sure why we do what appears to be a similar 1260 operation twice. It also seems odd that this is in the CombineFinal 1261 step, and not in the Reject step.} Finally, the rejected pixels are 1262 allowed to grow to include pixels that are neighbors to many rejected 1263 pixels. The ISIS kernel used in the previous step is used to 1264 determine the largest square box that contains under the limit of 1265 $0.25 * \sum_{x,y} kernel**2$. This box is then convolved with the 1266 rejected pixel mask to reject their neighbors. 1267 1268 This final list of rejected pixels is passed to the final combination 1269 pass, which does not iterate, and simply excludes the rejected 1270 pixels. \czwdraft{This is a bad paragraph.} 1271 1272 \czwdraft{We make the stacked image, the stacked variance, the stack 1273 mask, the exposure time mask, the exp weight containing the weighted 1274 exposure times, and a number image, containing the number of inputs 1275 used for each pixel.} 1276 1277 1278 1279 1280 1281 % FINAL COMBINE 1282 % Grow rejected pixels 1283 %% set threshold of (POOR.FRACTION = 0.25) * sum(kernel)**2 1284 %% Choose the largest square box that contains just under that threshold. 1285 %% Convolve that box with the rejected pixels to grow them. 1286 % Run combination pass again, but without doing rejection, simply applying the rejection lists already calculated. 1287 % :: 1288 % if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) { 1289 % iter = 0 1290 % combineRej = NAN 1291 % combineSys = NAN 1292 % combineDiscard = NAN 1293 % if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej, 1294 % combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) { 1295 1296 1297 The convolved stack products are not retained, as the convolution 1298 reduces the resolution of the final image. Instead, we apply the 1299 normalizations and rejected pixel maps generated from the convolved 1300 stack process to the original unconvolved input images. This produces 1301 an unconvolved stack that has the optimum image quality possible from 1302 the input images. Not convolving does mean that the PSF shape changes 1303 somewhat across the image, as the different FWHM of the input images 1304 print through in the different regions in which they have contributed 1305 to the final image. 1306 1307 % UNCONVOLVED IMAGE 1308 % if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false, true)) { 1309 % no grow 1310 1311 % only retain unconvolved products. 1312 1313 1314 One benefit of producing the final stacked image from the weighted 1315 mean of the unrejected input images is that faint sources do not have 1316 their contribution removed as much as median filtering would allow. 1317 \czwdraft{I did something to prove this once, but can't find it right 1318 now. Comparing the ppStack output catalog to one constructed from a 1319 simple median filtered stack shows that the ppStack catalog detects 1320 sources up to 0.XX magnitudes fainter than the median stack. This 1321 does increase the possibility of false positives.} 750 1322 751 1323 … … 755 1327 756 1328 \end{document} 1329 1330 1331 % http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation
Note:
See TracChangeset
for help on using the changeset viewer.
