Index: trunk/doc/release.2015/ps1.detrend/detrend.tex =================================================================== --- trunk/doc/release.2015/ps1.detrend/detrend.tex (revision 39231) +++ trunk/doc/release.2015/ps1.detrend/detrend.tex (revision 39232) @@ -139,22 +139,20 @@ \section{INTRODUCTION}\label{sec:intro} - +%% http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?2007ASPC..364..153M&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf \section{Camera description} -\czwdraft{reference to original paper} - -\czwdraft{60 otas} - -\czwdraft{64 cells per ota} - -\czwdraft{effectively 60x64 different cameras, each with particular gain/noise/etc characteristics} - -\czwdraft{Add summary of detrending steps} - -\czwdraft{Summary of detrending steps with references to the sections} +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. + +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}. + +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. + + + +% Discuss 2-phase/3-phase device differnces \section{Burntool / Persistence effect} -Stars that are nearing saturation \czwdraft{(30000 DN)} cause +Stars that are nearing saturation on GPC1 cause persistance problems during the read out of the image, creating trails of light are left on the image. During the read out process of an @@ -179,6 +177,5 @@ BURNTOOL program. This program does an initial scan of the images, and identifies stars brighter than a given threshold of 30000 DN. The -trail from that star is fit with a one-dimensional power law -\czwdraft{in each pixel column}, based on empirical evidence that this +trail from that star is fit with a one-dimensional power law in each pixel column, based on empirical evidence that this is the functional form of this persistence effect. Once this fit is done, the model is subtracted from the image, and the location of the @@ -218,10 +215,11 @@ \section{Masking} +\czwdraft{Technically, we mask the image prior to burntool application now.} \subsection{Static Masks} Due to the large size of the detector, it is to be expected that there -will be a number of pixel defects that \czwdraft{do not measure light} -as well as their neighbors. To remove these pixels, we have +will be a number of pixel defects that do not have the detection sensitivity on par +with their neighbors. To remove these pixels, we have constructed a static mask that identifies the known defects. This mask is constructed in three phases. @@ -231,31 +229,50 @@ detector. Twenty-five of the sixty OTAs in GPC1 show some evidence of CTE issues, with this pattern showing up (to varying degrees) in -triangular sets of cells on the OTA. \czwdraft{probably a figure would - help explain this?} To generate the mask, a sample set of evenly -illuminated flat field images are measured to produce a map of the -image variance in 20x20 pixel bins. As the flat image largely -illuminates the image uniformly, the expected variances should be -Poissonian distributed with the flux level. However, in regions with -CTE issues, adjacent pixels are not independent, allowing the charge -to spread. This reduces the pixel-to-pixel differences, resulting in -a lower-than-expected variance. All regions with variance -\czwdraft{X} smaller than expected are added to the static CTEMASK. - -The next step of mask construction is to examine the detector for +roughly triangular patches on the OTA due to defects in the +semiconductor \czwdraft{doping}. To generate the mask for these +regions, a sample set of \czwdraft{N} evenly illuminated flat field +images were measured to produce a map of the image variance in 20x20 +pixel bins. As the flat image is expected to illuminate the image +uniformly, the expected variances in each bin should be Poissonian +distributed with the flux level. However, in regions with CTE issues, +adjacent pixels are not independent, allowing the charge in those +pixels to spread. This reduces the pixel-to-pixel differences, +resulting in a lower-than-expected variance. All regions with +variance \czwdraft{0.5} smaller than expected are added to the static +CTEMASK. + +The next step of mask construction is to examine the flat and dark +models, and exclude pixels that appear to be poorly corrected by these +models. The darkmask process looks for pixels that are more than +\czwdraft{8} sigma discrepant in \czwdraft{10\%} of the +\czwdraft{test} images after those images have had the dark model +applied to them. These pixels are assumed to be unstable with respect +to the dark model, and have the DARK bit set in the static mask, +indicating that they are unreliable in scientific observing. +Similarly, the flatmask process looks for pixels that are \czwdraft{3} +sigma discrepant in the same fraction of \czwdraft{test} images after +both the dark and flat models have been applied. Those pixels that do +not follow the flat field model of the rest of image are assigned the +FLAT mask bit in the static mask, removing the pixels that cannot be +corrected to a linear response. + +The final step of mask construction is to examine the detector for bright columns and other static pixel issues. This is first done by -\czwdraft{I think Heather wrote a program to do this, but I'm not - totally sure how it works} scanning a set of images for pixels that -have values that do not change throughout a sequence of \czwdraft{N} -exposures. Such common pixel values cannot be caused by astronomical -effects, and must be due to the detector itself. This does an -excellent job of removing the majority of the problem pixels. A -manual inspection allows human interaction to identify other -inconsistent pixels including the vignetted regions around the edge of -the detector. \czwdraft{This might be a lie} As the size of the -vignetted region changes with filter, we have taken the g filter as -the baseline to define the static mask, resulting in the smallest -possible unvignetted region. - -The final static mask is the union of the CTE mask, the manual mask, \czwdraft{make this a paragraph}. +processing a set of \czwdraft{100 i filter} science images in the same +fashion as for the darktest. A median image is constructed from these +inputs along with the per-pixel variance. These images are used to +identify pixels that have unexpectedly low variation between all +inputs, as well as those that significantly deviate from the global +median value. Once this initial set of bad pixels is identified, a +$3\times{}3$ pixel triangular kernel is convolved with the initial +set, and any convolved pixel with value greater than \czwdraft{1.0} is +assigned to the static mask. This does an excellent job of removing +the majority of the problem pixels. A subsequent manual inspection +allows human interaction to identify other inconsistent pixels +including the vignetted regions around the edge of the detector. +\czwdraft{This might be a lie} As the size of the vignetted region +changes with filter, we have used the g filter to set the baseline +unvignetted region to define the static mask, resulting in the +smallest possible unvignetted region. \begin{figure} @@ -268,9 +285,11 @@ defects, we also generate a set of dynamic masks that change with the astronomical features in the image. These masks are advisory in -nature, and no not completely exclude the pixel from further -consideration. The first of these dynamic masks indicates the -presence of a corrected burntool trail. These pixels are included for -phtometry, but are rejected more readily in the stacking and -difference image construction. +nature, and do not completely exclude the pixel from further +processing consideration. The first of these dynamic masks indicates +the presence of a corrected burntool trail. These pixels are included +for phtometry, but are rejected more readily in the stacking and +difference image construction, as they are more likely to have small +residual contributions from the under or over subtraction of the +burntool correction. The remaining dynamic masks are not generated until the IPP camera @@ -278,23 +297,68 @@ photometry is complete, and an astrometric solution is known for the exposure. This added information provides the positions of bright -sources, which are the origin for the image artifacts that the dynamic -mask identifies. +sources based on the reference catalog, including those that fall +slightly out of the detector field of view or within the inter chip +gaps, where internal photometry may not have identified them. These +bright sources are the origin for many of the image artifacts that the +dynamic mask identifies and excludes. + + +\begin{deluxetable}{ccl} + \tablecolumns{3} + \tablewidth{0pc} + \tablecaption{GPC1 Mask Values} + \tablehead{\colhead{Mask Name} & \colhead{Mask Value} & \colhead{Description}} + \startdata + DETECTOR & 0x0001 & A detector defect is present. \\ + FLAT & 0x0002 & The flat field model does not calibrate the pixel reliably. \\ + DARK & 0x0004 & The dark model does not calibrate the pixel reliably. \\ + BLANK & 0x0008 & The pixel does not contain valid data. \\ + CTE & 0x0010 & The pixel has poor charge transfer efficiency. \\ + SAT & 0x0020 & The pixel is saturated. \\ + LOW & 0x0040 & The pixel has a lower value than expected. \\ + SUSPECT & 0x0080 & The pixel is suspected of being bad. \\ + BURNTOOL & 0x0080 & The pixel may contain an uncorrected or over-corrected burntool streak. \\ + CR & 0x0100 & A cosmic ray is present. \\ + SPIKE & 0x0200 & A diffraction spike is present. \\ + GHOST & 0x0400 & An optical ghost is present. \\ + STREAK & 0x0800 & A streak is present. \\ + STARCORE & 0x1000 & A bright star core is present. \\ + CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\ + CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\ + MARK & 0x8000 & An internal flag for temporarily marking a pixel. \\ + \enddata + \label{tab:mask_values} +\end{deluxetable} + \subsubsection{Crosstalk ghosts} + Due to electrical crosstalk between the flex cables connecting the -individual detectors, ghost objects can be created on some OTAs due to -the presence of a bright object in a different position. Table -\ref{tab:crosstalk_rules} summarizes the list of known crosstalk -rules. In each of these cases, a source object brighter than -14.47 -magnitude (instrumental) creates a ghost object many orders of -magnitude fainter at the target location. The cell (x,y) coordinate -is identical between source and ghost, as a result of the transfer -occurring as the devices are read. A circular mask is asdded to the -ghost location with radius $R = 3.44 \left(-14.47 - m_{source, - instrumental}\right)$. Any objects in the photometric catalog found -at the location of the ghost mask have a flag set, marking the object -as a ghost. - -\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.} +individual detector devices, ghost objects can be created on some OTAs +due to the presence of a bright source at a different position on the +camera. Table \ref{tab:crosstalk_rules} summarizes the list of known +crosstalk rules. In each of these cases, a source object brighter +than -14.47 magnitude (instrumental) creates a ghost object many +orders of magnitude fainter at the target location. The cell (x,y) +coordinate is identical between source and ghost, as a result of the +transfer occurring as the devices are read. A circular mask is asdded +to the ghost location with radius $R = 3.44 \left(-14.47 - m_{source, + instrumental}\right)$ pixels. Any objects in the photometric +catalog found at the location of the ghost mask have a \czwdraft{flag} +set, marking the object as a likely ghost. The majority of the +crosstalk rules are bi-directional, with a source in either position +creating a ghost at the corresponding crosstalk target position. The +two faintest rules are uni-directional, likely due to differences in +the \czwdraft{magical properties of the electronics}. + +For the very brightest sources ($m_{instrumental} < -15$), there can +be crosstalk ghosts between all columns of cells during the readout. +These ``bleed'' ghosts were originally identified as ghosts of the +saturation bleeds appearing in the neighboring cells, and as such, the +masking for these objects puts a rectangular mask down from top to +bottom of cells in all columns that are in the same row of cells as +the bright source. The width of this box is a function of the source +magnitude, with $W = 5 * \left(-15 - m_{source, instrumental}\right)$ + pixels. \begin{deluxetable}{lllc} @@ -318,23 +382,27 @@ \end{deluxetable} +\begin{figure} + \caption{Figure of crosstalk ghost and bright star source. Plot of cut across ghost to illustrate the flat-top shape.} +\end{figure} \subsubsection{Optical ghosts} - -Due to an issue with the anti-reflective coating, bright sources can -also result in large out of focus objects, particularly in the +% http://arxiv.org/pdf/1207.2513v1.pdf +Due to imperfections in the anti-reflective coating, bright sources +can also result in large out of focus objects, particularly in the g-filter data. These objects are the result of light reflecting back -off the surface of the detector, reflecting again off the \czwdraft{No - clue} mirror, and then back down onto the focal plane. Due to the -extra travel distance, the resulting source is out of focus and -elongated along the radial direction of the telescope. These optical -ghosts can be modeled as a bright star in location (X,Y) on the focal -plane creates a reflection ghost on the opposite side of the optical -axis at (-X,-Y). The exact location is fit as a third order -polynomial in the focal plane x and y directions. An elliptical -annulus mask is constructed at the expected ghost location, with the -major and minor axes defined by linear functions of the ghost distance -from the optical axis, and orientation \czwdraft{pointing along - radius}. All stars brighter than a filter-dependent threshold -(listed in table \ref{tab:ghost_magnitudes}) have masks constructed. +off the surface of the detector, reflecting again off the lower +surfaces of the optics (particularly the L1 corrector lens), and then +back down onto the focal plane. Due to the extra travel distance, the +resulting source is out of focus and elongated along the radial +direction of the telescope. These optical ghosts can be modeled as a +bright star in location (X,Y) on the focal plane creates a reflection +ghost on the opposite side of the optical axis at (-X,-Y). The exact +location is fit as a third order polynomial in the focal plane x and y +directions. An elliptical annulus mask is constructed at the expected +ghost location, with the major and minor axes defined by linear +functions of the ghost distance from the optical axis, and oriented +along the radius of the detector. All stars brighter than a +filter-dependent threshold (listed in table +\ref{tab:ghost_magnitudes}) have such masks constructed. \begin{deluxetable}{lc} @@ -356,65 +424,39 @@ \czwdraft{include full polynomial forms? How best to do that?} - -%% -%% GHOST.CENTER.X METADATA -%% NORDER_X S32 3 -%% NORDER_Y S32 3 -%% VAL_X00_Y00 F64 -1.215661e+02 -%% VAL_X01_Y00 F64 1.321875e-02 -%% VAL_X02_Y00 F64 -4.017026e-09 -%% VAL_X03_Y00 F64 1.148288e-10 -%% VAL_X00_Y01 F64 -1.908074e-03 -%% VAL_X01_Y01 F64 8.479150e-08 -%% VAL_X02_Y01 F64 1.635732e-11 -%% VAL_X00_Y02 F64 2.625405e-08 -%% VAL_X01_Y02 F64 1.125586e-10 -%% VAL_X00_Y03 F64 2.912432e-12 -%% NELEMENTS S32 10 -%% END - -%% GHOST.CENTER.Y METADATA -%% NORDER_X S32 3 -%% NORDER_Y S32 3 -%% VAL_X00_Y00 F64 2.422174e+01 -%% VAL_X01_Y00 F64 4.170486e-04 -%% VAL_X02_Y00 F64 -1.934260e-08 -%% VAL_X03_Y00 F64 -1.173657e-12 -%% VAL_X00_Y01 F64 1.189352e-02 -%% VAL_X01_Y01 F64 -9.256748e-08 -%% VAL_X02_Y01 F64 1.140772e-10 -%% VAL_X00_Y02 F64 8.123932e-08 -%% VAL_X01_Y02 F64 1.328378e-11 -%% VAL_X00_Y03 F64 1.170865e-10 -%% NELEMENTS S32 10 -%% END -%% # These are the original linear solutions -%% GHOST.INNER.MAJOR METADATA -%% NORDER_X S32 1 -%% VAL_X00 F64 3.926693e+01 -%% VAL_X01 F64 5.325759e-03 -%% NELEMENTS S32 2 -%% END - -%% GHOST.INNER.MINOR METADATA -%% NORDER_X S32 1 -%% VAL_X00 F64 5.287548e+01 -%% VAL_X01 F64 -2.191669e-03 -%% NELEMENTS S32 2 -%% END - -%% GHOST.OUTER.MAJOR METADATA -%% NORDER_X S32 1 -%% VAL_X00 F64 7.928722e+01 -%% VAL_X01 F64 1.722181e-02 -%% NELEMENTS S32 2 -%% END - -%% GHOST.OUTER.MINOR METADATA -%% NORDER_X S32 1 -%% VAL_X00 F64 1.314265e+02 -%% VAL_X01 F64 -2.627153e-03 -%% NELEMENTS S32 2 -%% END +\begin{deluxetable}{lcc} + \tablecolumns{3} + \tablewidth{0pc} + \tablecaption{Optical Ghost Center Transformations} + \tablehead{\colhead{Polynomial Term}&\colhead{X center}&\colhead{Y center}} + \startdata + $x^0 y^0$ & -1.215661e+02 & 2.422174e+01 \\ + $x^1 y^0$ & 1.321875e-02 & 4.170486e-04 \\ + $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\ + $x^3 y^0$ & 1.148288e-10 & -1.173657e-12 \\ + $x^0 y^1$ & -1.908074e-03 & 1.189352e-02 \\ + $x^1 y^1$ & 8.479150e-08 & -9.256748e-08 \\ + $x^2 y^1$ & 1.635732e-11 & 1.140772e-10 \\ + $x^0 y^2$ & 2.625405e-08 & 8.123932e-08 \\ + $x^1 y^2$ & 1.125586e-10 & 1.328378e-11 \\ + $x^0 y^3$ & 2.912432e-12 & 1.170865e-10 \\ + \enddata + \label{tab:ghost_centers} +\end{deluxetable} + +\begin{deluxetable}{lcccc} + \tablecolumns{5} + \tablewidth{0pc} + \tablecaption{Optical Ghost Annulus Axis Length} + \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}& \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}} + \startdata + $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\ + $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\ + \enddata + \label{tab:ghost_radii} +\end{deluxetable} + +\begin{figure} + \caption{Figure of full FOV showing optical ghosts. Possibly only a few OTAs to illustrate shape deformation.} +\end{figure} \subsubsection{Glints} @@ -429,6 +471,6 @@ have a dynamic mask constructed when a reference source falls on the focal plane within \czwdraft{approximately} one degree of the detector -edge. This mask is 150 pixels wide, and $L = 2500 \left(-20 - -m_{inst}\right)$. +edge. This mask is 150 pixels wide, with length $L = 2500 \left(-20 - +m_{inst}\right)$. \czwdraft{Am I correct that this is basically a one-degree edge around the detector?} %% @@ -456,4 +498,8 @@ %% END +\begin{figure} + \caption{Example of glint.} +\end{figure} + \subsubsection{Diffraction spikes} @@ -462,5 +508,5 @@ with length $L = 10^{0.096 * (7.35 - m)} - 200$ and width $W = 8 + (L - 200) * 0.01$. These spikes are dependent on the camera rotation, -and are oriented at $\theta = n * \frac{pi}{2} - \mathrm{ROTANGLE} + +and are oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} + 0.798$. @@ -469,5 +515,7 @@ The cores of saturated stars are masked as well, with radius $r = 10.15 * (-15 - m_{inst})$. \czwdraft{good job here.} -\czwdraft{Write up something about the masking fraction.} +\begin{figure} + \caption{Example of saturated star, which will also nicely show the diffraction spikes.} +\end{figure} \subsection{Video Mask} @@ -493,13 +541,27 @@ \subsection{Masking fraction} -\czwdraft{\% due to chip/cell gaps} - -\czwdraft{\% due to faulty pixels} - -\czwdraft{\% due to CTE} - -\czwdraft{\% due to vinetting} - -\czwdraft{\% average dynamic masking} +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\%. + +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. + +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\%. + +%% 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; +%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+ +%% | 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) | +%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+ +%% static dynamic advisory +%% | g.00000 | 0.19642137972007 | 0.00010322263512709 | 0.026838445469766 +%% | 0.20949461794863 | 9.89200027293e-05 | 0.026431927734548 | +%% | r.00000 | 0.19675996201399 | 0.00025214447869606 | 0.032641054600788 +%% | 0.20989768279138 | 0.00023994155711801 | 0.032178525485201 | +%% | i.00000 | 0.19677587604327 | 0.00057470697316504 | 0.038096251937072 +%% | 0.21003570722292 | 0.00053987093278142 | 0.037471018638997 | +%% | z.00000 | 0.1974290315691 | 0.00024758901226967 | 0.03064123748973 +%% | 0.21055007930696 | 0.00023452690039757 | 0.030144453360769 | +%% | y.00000 | 0.19828990634315 | 0.00014523787521897 | 0.021984846417987 +%% | 0.21130344126869 | 0.00013634812877977 | 0.02163070300815 | + +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. \section{Overscan} @@ -524,6 +586,6 @@ frames with a ramp of exposure times. As the exposure time increases, the flux on each pixel also increases in what is expected to be a -linear manner. Each of these dark ramp exposures is overscan -corrected, and then the median is calculated for each cell, as well as +linear manner. Each of these dark exposures in this exposure time ramp is overscan +corrected, and then the median is calculated for each cell, as well as for the rows and columns within ten pixels of the edge of the science region. From these median values at each exposure time value, we can @@ -537,22 +599,32 @@ cell results in the center of the detector. -This non-linearity effect appears to be stable in time, with no -evident change over a year's worth of data. +This non-linearity effect appears to be stable in time, with little +evident change over the survey duration. \czwdraft{I have figures at http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity that might be useful} +%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity_AllEdges +%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearityArchive + +\begin{figure} + \caption{Example plot of linearity as a function of incident brightness.} +\end{figure} \section{Dark/Bias Subtraction} - +% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Background_Dark_Model The dark model we make for GPC1 considers each pixel individually, -independent of any neighbors. To create the dark model for each -pixel, we fit an arbitrary dimensional model \czwdraft{clunky} to the -array of input pixels from a selection of dark images. The current -model is linear \czwdraft{really?} in both the exposure time and the -detector temperature. Adding in a constant value for the fit provides -three parameters that define the dark model for that pixel. As this -constant value is effectively the bias value for that pixel, we do not -do a separate bias correction. This model is applied to science -images by fitting the correct dark value based on the exposure time -and detector temperature for that exposure. +independent of any neighbors. To create the dark model, we fit an multi-dimensional model to the array of input pixels +from a randomly selected set of 100-150 \czwdraft{overscan corrected} +dark frames chosen from a given date range. The model fits +each pixel as a function of the exposure time $t_{exposure}$ and the +detector temperature $T_{chip}$ such that $dark = a_0 + a_1 +t_{exposure} + a_2 T_{chip} t_{exposure} + a_3 T_{chip}^2 +t_{exposure}$. This fitting is performed over the sample of input pixels, +and the coefficients $a_i$ stored in the detrend image. The constant +$a_0$ term includes the bias signal, and as such, a separate bias +subtraction is not necessary. + +Applying the dark model is simply a matter of calculating the response +for the exposure time and detector temperature for the image to be +corrected, and subtracting the resulting dark signal from the image. \subsection{Time evolution} @@ -560,8 +632,7 @@ \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.} -Unfortunately, the dark model is not consistently stable on the time -span of multiple months. Some of the changes in the dark can be +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 attributed to changes in the voltage settings of GPC1, but the -majority seem to be the result of some unknown parameter. Largely, we +majority seem to be the result of some unknown parameter. We can separate the dark model history of GPC1 into three epochs. The first epoch covers all data taken prior to 2010-01-23. This epoch @@ -572,5 +643,5 @@ characterized by a largely stable but oscillatory dark solution. There appear to be two modes that the dark model switches between -apparently at random. No clear cause has been established for these +apparently at random. No clear cause has been established for the switching, but there are clear differences between the two modes \czwdraft{figures?}. @@ -586,5 +657,7 @@ appropriate ``A'' and ``B'' mode dark frames. Using the appropriate dark minimizes the effect of this bias gradient in the dark corrected -data. +data. Table \ref{tab:dark mode dates} lists the dates used for each dark mode. + +\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.} After 2011-05-01, the two-mode behavior of the dark disappears, and is @@ -594,16 +667,26 @@ These darks cover the range from 2011-05-01 to 2011-08-01, 2011-08-01 to 2011-11-01, and 2011-11-01 and on. The reason for this time -dependent drift is unknown, but we seem to be able to model it with +evolution is unknown, but we seem to be able to model it with reasonable accuracy by creating new dark models. +\begin{figure} + \caption{Example of raw and dark calibrated exposure. Plots of horizontal cuts for A/B/average corrections.} +\end{figure} + +\subsection{Video Dark} + +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. + +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}$. + \section{Noisemap} -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 read out. 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. - -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 teh 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. - -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}$. - -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. +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. + +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. + +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}$. + +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. \section{Remnance?} @@ -623,5 +706,5 @@ uniformly illuminated image. Using a dome screen is not possible, as the variations in illumination and screen rigidity create unusably -large scatter between different images. Because of this, we use sky +large scatter between different images that are caused by the detector response function. Because of this, we use sky flat images taken at twilight, which are more consistently illuminated than screen flats. We calculate the mean of these images to determine @@ -630,12 +713,15 @@ From this initial flat model, we construct a correction to remove the effect of the problems illuminating the large area. This is done by -dithering a series of exposures across a given pointing. By comparing -the measured fluxes for a given star as a function of position, we can -correct out the errors in the flat model. +dithering a series of science exposures across a given pointing. By +comparing the measured fluxes for a given star as a function of +position on the detector, we can determine the position dependent +scaling factors. These scale factors can then be used to correct the +initial flat field model to better represent the detector response. The flat model appears stable with time, although directly measuring this is as difficult as originally constructing the model. However, -due to the photometric consistency observed in GPC1 measurements, we -can be confident that the flat model is not changing much. +due to the photometric consistency observed in the catalog of GPC1 measurements, we +can be confident that the flat model is not as time dependent as the +dark correction. @@ -647,5 +733,5 @@ to remove the remaining row-by-row variation, and the PATTERN.CELL and PATTERN.CONTINUITY corrections attempt to ensure that the cells of a -given OTA are consistent with each other. These corrections are +given OTA are consistent with the other cells on that OTA. These corrections are largely designed to fix issues that are not stable enough with time for the dark model or flat field model to fully account for the @@ -653,15 +739,19 @@ \subsection{Pattern Row} - +% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study As discussed above in the dark and noisemap sections, certain -detectors have significant row-by-row bias offsets. As the level of -the offset is largely random, the dark correction cannot fully remove -this structure from the images. Therefore, we apply the PATTERN.ROW -correction in an attempt to mitigate the offsets. To force the rows -to agree, a \czwdraft{first} order polynomial is fit to each row in the -cell, and that trend subtracted from the data. The median offset -(corresponding to the background level) is then added back to the -image so that the cell matches its neighbors during background -subtraction. +detectors have significant row-by-row bias offsets. The magnitude of +these offsets increases as the distance from the readout amplifier +increases, resulting in horizontal streaks that are more pronounced +along one edge of the cell. As the level of the offset is largely +random, the dark correction cannot fully remove this structure from +the images, and the noisemap level only indicates the level of the +variance added by these bias offsets. Therefore, we apply the +PATTERN.ROW correction in an attempt to mitigate the offsets. To +force the rows to agree, a \czwdraft{second} order polynomial is fit to +each row in the cell, and that trend subtracted from the data. The +median offset (corresponding to the background level) for each row is then fit by a first order polynomial, and that trend is then added +back to the image so that the sky level on the cell matches its neighbors during +background subtraction. This correction was required on all cells on all OTAs prior to @@ -669,5 +759,5 @@ electronics resolved the row-by-row offsets for the majority of the detectors. As a result, we only apply this correction where it is -necessary, as shown in figure \czwdraft{X}. +necessary, as shown in Figure \ref{fig: pattern row required}. Although this correction does resolve the row-by-row offset issue in a @@ -676,15 +766,19 @@ near these objects. As the offsets are calculated on the pixel rows, this oversubtraction is not uniform around the object, but is -preferentially along the $\pm x$ axis of the object. - -\czwdraft{keep this?} This row-by-row offset is visible in similar -camera designs, and has been removed by identifying the noise signal -in the pixel data stream. By taking the FFT of the pixels and a -reference signal, the frequency of this noise can be isolated and -removed, resulting in a much cleaner image. However, GPC1 does not -record the value of the reference signal, instead automatically -subtracting it from the data values. Without this comparison signal, -we have been unable to reproduce this method, as there is no obvious -FFT component visible. +preferentially along the horizontal x axis of the object. + +%% \czwdraft{keep this?} This row-by-row offset is visible in similar +%% camera designs, and has been removed by identifying the noise signal +%% in the pixel data stream. By taking the FFT of the pixels and a +%% reference signal, the frequency of this noise can be isolated and +%% removed, resulting in a much cleaner image. However, GPC1 does not +%% record the value of the reference signal, instead automatically +%% subtracting it from the data values. Without this comparison signal, +%% we have been unable to reproduce this method, as there is no obvious +%% FFT component visible. + +\begin{figure} + \caption{Example of pre/post pattern row application.} +\end{figure} \subsection{Pattern Cell} @@ -692,37 +786,37 @@ As the bias level of a given cell may not exactly match that of its neighbors, fitting a smooth background model results in over and -under-subtraction of the sky level at these discontinuities. The -PATTERN.CELL correction was the first attempt to remove this effect on -the worst cells, by forcing all the cells of an OTA to the same level. -Each cell has the median value measured, and then an offset added that -shifts each cell to match the median of these medians. +under-subtraction of the sky level at the cell boundary +discontinuities. The PATTERN.CELL correction was the first attempt to +remove this effect on the worst cells, by forcing all the cells of an +OTA to the same level. Each cell has the median value measured, and +then each cell has an offset added that shifts the cell to match the +median of those medians. This correction is reasonable when the astronomical signal is smooth, with no objects that are large relative to the size of an individual cell. However, the presence of large galaxies (or even bright stars) -can force some cells into a nearly arbitrary offset from their -neighbors. Because of this issue, we no longer apply this correction -to any data. +can bias the offsets for some cells from their neighbors. Because of +this issue, we no longer apply this correction to any data. \subsection{Pattern Continuity} -As the PATTERN.CELL correction was clearly defective in many -situations, we designed a replacement correction that would distort -large objects less. In addition, studies of the background level -illustrated that the row-by-row bias introduces a small background -gradient along the rows of the cells. This results in a ``sawtooth'' -pattern across an OTA, and as the background model assumes a smooth -sky level, we saw evidence of over and under subtraction at cell -boundaries. As the PATTERN.CELL was designed to correct mean changes -between cells, it could not adequately resolve this higher order -issue. +As the PATTERN.CELL correction was clearly insufficient in many +situations, we designed a replacement correction that would lower the +distortion for large objects less. In addition, studies of the +background level illustrated that the row-by-row bias introduces +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 +in a ``sawtooth'' pattern horizontally across an OTA, and as the background model +assumes a smooth sky level, this induces over and under +subtraction at cell boundaries. As the PATTERN.CELL was designed to +correct mean changes between cells, it could not adequately resolve +this higher order issue. The replacment for PATTERN.CELL was the PATTERN.CONTINUITY correction, which attempts to match the edges of a cell to those of its neighbors. -For each cell, a thin box on each edge is extracted and the median +For each cell, a thin box \czwdraft{10} pixels wide on each edge is extracted and the median value calculated for that box. These median values are then used to construct a vector of differences $diff_i = \sum_{j,j'} Edge_{i,j} - Edge)_{i',j'}$, along with a matrix of associations $A_{i,i'} = -\sum_{j,j'} \delta(j,j')$ denoting which cell boundary touches +\sum_{j,j'} \delta(j,j')$ denoting which cell boundaries touch another. By solving the system $A x = diff$, we can find the set of offsets $x_i$ that should be applied to each cell to ensure the @@ -732,20 +826,498 @@ to align the cells into a single ramp, at the expense of the absolute background level. However, as we subtract off a smooth background -model, this absolute level is unimportant. The fact that the final +model, the deviations from an absolute sky level are unimportant. The fact that the final ramp is smoother than it would be otherwise also allows for the background subtracted image to more closely match the astronomical sky, without over- and under-subtractions at cell edges. -%% \section{Fringe correction} - -%% \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.} - -%% \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.} - -%% \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.} - -%% \section{Background subtraction} - -%% \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.} +\begin{figure} + \caption{Continuity example, with background issue.} +\end{figure} + +\section{Fringe correction} + +Due to variations in the thickness of the detectors, we observe +interference patterns at the infrared end of the filter set, as +the wavelength of the light becomes comparable to the thickness of these variations. +Visually inspecting the images shows that the fringing is most +prevalent in the y-filter images, with minimal fringing in other +bands. As a result of this, we only apply a fringe correction to the y filter +data. + +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 +filtered out as an additional noise term. A least squares fit to the +data is then calculated, providing the model of the fringe strength at +that location. + +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. + +\section{Background subtraction} + +\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.} + +%% * Magic +%% * Warping +%% * warping kernel +%% * linear-by-pieces +%% * Covariance +%% * def of skycells? +%% * Stacking +%% * pixel combination rules +%% * pixel rejections +%% * convolution for matching (success and failure) +%% * Difference Image analysis + + +\section{Warping} + +To provide a consistent and uniform set of images for co-added image +stacking and image differences, the individual mosaicked OTA images +are projected onto a common set of tangent plane projected regions. +These projection cells are $4\times{}4$ degree fields spaced onto set +of projection centers that fully cover the sky. These projection +cells are arranged into rings of constant declination, and allowed to +overlap as $|\delta|$ increases. Each projection cell is further +subdivided into \czwdraft{size} sky cells, which have constant overlap +regions of \czwdraft{overlap}. These skycells are the main image unit +used for processing image data beyond the initial chip stage. + +After the detrending and photometry, the detection catalog for the +full camera is fit to the reference catalog, producing third-order +astrometric solutions that map the detector focal plane to the sky, +and map the individual OTA pixels to the detector focal plane. This +solution is then used to determine which skycells the exposure OTAs +overlap. + +Foreach output skycell, all overlapping OTAs and the calibrated +catalog are read into the \textbf{pswarp} program. Each input image +is examined in order, and the same transformation performed. This +transformation breaks the output warp image into $128\times{}128$ +pixel grid boxes. Each grid box has a locally linear map calculated +that converts the output warp image coordinates to the input chip +image coordinates. By doing the transformation in this direction, +each output pixel has a unique sampling position on the input image +(although it may be off the image frame and therefore not populated), +preventing gaps in the output image due to the spacing of the input +pixels. + +With the locally linear grid defined, Lanczos interpolation with +filter size parameter $a = 3$ on the input image is used to determine +the values to assign to the output pixel location. The output +locations are shifted by 0.5 pixels to let the interpolation select +the value that would be assigned to the center of the output +pixel. This process is repeated for all grid boxes, for all input +images, and for each output image product: the science image, the +variance, and the mask. \czwdraft{The jacobian is multiplied to the + image value, and squared and multiplied to the variance. I don't + understand that.} + +As the interpolation constructs the output pixels from more than one +input pixel, there is a covariance term that is must be included. For +each locally linear grid box, the covariance is calculated from the +kernel in the center of the 128 pixel range. Once the image has been +fully populated, this set of individual covariance matrices is +averaged to create the final covariance for the full image. + +An output catalog is also constructed from the full exposure input +catalog, including only those objects that fall on the warped image. +These detections are transformed to match the new image location, and +to scale the position errors based on the new orientation. + +The output image also contains header keywords SRC\_0000, SEC\_0000, +MPX\_0000, and MPY\_0000 that contain the mappings from the warped +pixel space to the input image. The SRC keyword lists the input OTA +name, and the SEC keyword lists the image section corresponding to the +locally linear grid box. The MPX and MPY contain the transformation +parameters for the locally linear grid. \czwdraft{Is this accurate?} + +% Read all images and astrometry +% Check which input images overlap with output image. => 8007 when the inputs don't overlap. +% Loop over each image. +% Append detections from input into output detection list. +% Determine transform back from warp pixels to source image. +%% 2nd order polynomial in both x and y for this transformation. and save to header +% Break warp image into 128x128pixel locally linear areas +% Mask non finite pixels as saturated. +% Define Lanczos-3 interpolation over the input image. +% Iterate over warp pixel space (on each locally linear grid) and map interpolated input pixel positions onto warp. +% Warp pixel space is defined as center based, so that's where the intpolation point comes from. +% Covariance calculated based on the interpolation kernel at the center of the ll grid. +% image = interp_image * jacobian +% var = interp_var * jacobian**2 +% mask = interp_mask +% jacobian = abs(mapXx * mapYy - mapYx * mapXy) +% I don't understand that jacobian. +% + + +\section{Stacking} + +Once individual exposures have been warped onto a common projection +system, they can then be combined without that added concern. In +order to obtain detections of faint images, and to provide a static +sky image without transient features, we coadd the individual warps +into a stacked image. Creating this stack also allows a complete +image to be constructed that does not have regions masked due to +falling between devices. + +The stacked image is comprised of all warp frames for a given skycell +in a single filter. The source catalogs and images are loaded into +the \textbf{ppStack} program to do prepare the inputs and stack the +frames while rejecting bad pixels. + +Once all files are ingested, the first step is to measure the size and +shapes of the input image PSFs. We exclude images that have a PSF +FWHM greater than 10 pixels, as those images have the worst seeing and +would degrade the final output stack. A target PSF for the stack is +constructed from the envelope of all input PSFs, which sets the target +PSF at the largest value among the input PSFs for all radii. This PSF +is then circularized to prevent any of the input images from being +deconvolved when matched to the target. + +The input images also need to be normalized to prevent differences in +seeing and sky transparency from causing discrepancies during pixel +rejection. From the calibrated input catalogs, we have the +instrumental magnitudes of all sources, along with the airmass, image +exposure time, and zeropoint. All output stacks are calibrated to a +zeropoint of 25.0 in all filters, and to have an airmass of 1.0. The +output exposure time is set as the sum of the input exposure times. +With this information, we can determine the relative transparency for +each input image by comparing matched sources between the different +images. Each image then has a normalization factor defined, equal to +$norm_{image} = (ZP_{image} - ZP_{target}) - transparency_{image} - +2.5 * \log_{10} (t_{target} / t_{image}) - airmassTerm * +(airmass_{image} - airmass_{target})$. The input source catalog is +adjusted to reflect this normalization, which is also retained for +application when the pixels are combined. + +% PREPARE +% load sources +% load psfs +% determine target as envelope of input psfs +% FWHM clipping at 10 +% measure seeing +% - // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency +% // EAM : the discussion here was not quite right (or at least sloppy). Here is a replacement explanation: +% // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by: +% // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency +% // NOTE the sign of 'transparency' this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where +% // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app +% // we want to adjust the input images to be in a consistent flux system so that the +% // final stack can be generated with a specific target zero point. Any adjustment to +% // the flux scale of the image must be made in coordination with the resulting +% // zeropoint, exposure time, and airmass such that the above relationship yields the +% // same apparent magnitude for a given star: +% // m_inst_i : instrumental mags on input image (in) +% // m_inst_o : instrumental mags on re-normalized image (out) +% // 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 +% // 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 +% // zp_i, airmass_i, t_i, trans_i : reported or measured for input image +% // zp_o = zpTarget (from recipe) +% // airmass_o = airmassTarget (from recipe) +% // t_o = sumExpTime [sum of input exposure times: once images are scale to this time, they can be avereaged] +% // trans_o = 0.0 [obviously!] +% // we have 2 cases: (a) all reported ZPs are good or (b) some are bad: +% // (a) FPA.ZP = zp_i + c1 * airmass_i +% // --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i) +% // (b) zp[i] = c1 * airmass_i + 2.5log(exptime_i) +% // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong. +% // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i + +With the normalization factors and target PSF chosen, the convolution +kernels can be calculated for each image. ISIS kernels are used with +FWHM values of 1.5, 3.0, and 6.0 pixels and polynomial orders of 6, 4, +and 2. \czwdraft{Skipping this bit because I'm not completely sure I + understand it.} The image is then scaled by the normalization as +$renorm = 10^{-0.4 * norm_{image}} / norm_{convolution}$, and the + variance by the square of that value. + +% MATCH +% match to target PSF. +% use ISIS kernels to do matching/convolution +% Input sources used for psf matching. +% @ISIS.WIDTHS F32 1.5 3.0 6.0 # Gaussian kernel FWHM values +% @ISIS.ORDERS S32 6 4 2 # Polynomial orders for ISIS kernels + +Once the convolution kernels are defind for each image, they are used +to convolve the image to match the target PSF. Any input image that +has a $\chi^2$ value larger than 4.0$\sigma$ larger than the median +value is rejected from the stack. + +% CONVOLVE +% Normalization to match target zeropoint/exptime +% Reject images with bad match chi^2 values. MATCH.REJ * rms + median threshold. +% Additional variance from the convolution chi^2 +% Calculate image weights based on variance: W_i = 1 / (ROBUST_MEDIAN(variance_image_i) * CovarianceFactor) +% CovarianceFactor = covariance->kernel[0][0] + +Following the convolution, and initial stack is constructed. For a +given pixel coordinate, the values at that coordinate are extracted +from all input images. Images that have a suspect mask bit (including +the SUSPECT, BURNTOOL, SPIKE, STREAK, STARCORE, and CONV.POOR bit +values) are appended to a suspect pixel list for preferential +exclusion. Following this, the pixel values are combined and tested +to attempt to identify discrepant values that should be excluded. + +If only a single input is available, the initial stack contains the +value from that single input. If there are only two inputs, the +average of the two is used. These cases should occur only rarely in +the $3\Pi$ survey, as there are many input exposures that overlap any +particular point on the sky. The more common case for three or more +inputs constructs a weighted average from the inputs, with the weight +set as a single value for each input image, and defined as the inverse +of the median variance value from that image's associated variance +map. This weight is used for the image and the exposure weighted +image: + +\begin{eqnarray} + S_{value} &=& \sum_i\left(value_{i} * weight_i\right) / \sum_i\left(weight_i\right) \\ + S_{exp weight} &=& \sum_i \left(exptime_i * weight_i\right) / \sum_i\left(weight_i\right) \\ +\end{eqnarray} + +The pixel exposure time is simply the sum of the input exposure time values, and the output variance is + +\begin{eqnarray} + S_{variance} &=& 1 / \sum_i \left( 1 / variance_i \right) +\end{eqnarray} + +The output mask value is taken to be zero (no masked bits), unless +there were no valid inputs, in which case the BLANK mask bit is set. + +% INITIAL COMBINE +% Calculate weighted mean of input images +% mu = sum_i(f_i * W_i) / sum_i(W_i) +% sigma = 1 / sum_i(1 / W_i) +% nu = sum_i(m_i) +% // We're not using the input pixel variances to generate a weighted average for the pixel flux (because +% // that introduces systematic biases), so the variance of the output pixel value should simply be: +% // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 +% // This reduces, when the weights are all identically unity, to: +% // variance_combination = sum(variance_i) / N^2 +% // and if the variances are all equal: +% // variance_combination = variance_individual / N +% // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). +% sumValueWeight = sum_i(values * weights) +% sumWeight = sum_i(weights) +% sumVarianceWeight == sum( 1 / variances) +% sumExp = sum_i(exptimes) +% sumExpWeight = sum_i(exptims * weights) +% mean = sumValueWeight / sumWeight +% var = 1 / sumVarianceWeight +% exp = sumExp +% expWeight = sumExpWeight + +% EXCEPT: if N = 1, accept it. if N = 2, take average. + +% if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter, +% combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) { +%bool pmStackCombine( +% pmReadout *combined, // output stacked readout +% pmReadout *expmaps, // output exposure map information +% psArray *input, // input exposures +% psImageMaskType badMaskBits, // treat these bits as 'bad' +% psImageMaskType suspectMaskBits, // treat these bits as 'suspect' +% psImageMaskType blankMaskBits, // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?) +% int kernelSize, +% float iter, 0.5 +% float rej, 4.0 +% float sys, 0.1 +% float olympic, 0.2 +% bool useVariance, +% bool safe, +% int nminpix, +% bool rejection) +%{ + +% combineExtract +%% pixels with mask values as suspect are appended to suspect pixel list. +% combinePixels +%% As described above. + +Following this initial combination, a ``testing'' loop iterates in an +attempt to identify outlier points. Again, if only one input is +available, that input is accepted. If there are two inputs, $A$ and +$B$, then a check is made to see if $(0.5 * (value_A - value_B))**2 > +rej**2 * (variance_A + variance_B + (sys * value_A)**2 + (sys * +value_B)**2)$, where $rej$ is the number of sigma deviant a point +needs to be to be excluded, set to 4.0 for the PV3 processing, and +$sys$ is an estimate of the systematic error, taken to be 0.1. + + +\czwdraft{This discussion seems out of place, but I'm not sure where a + better place is.} Due to the various non-astronomical ghosts that +can occur on GPC1, and the fact that they may not be masked +aggressively enough to ensure all bad pixels are removed, it is +expected that some of the inputs for a given stack pixel are not in +agreement with the others. In general, there is the population of +input pixel values around the correct astronomical level, as well as +possible populations at lower pixel value (such as due to an +over-subtracted burntool trail) and at higher pixel values (such as +that caused by an incompletely masked optical ghost). Due to the +observation strategy to image a given field twice to allow for +warp-warp difference images to be constructed to identify transient +detections, higher pixel values that come from sources like optical +ghosts that are a function of pointing will come in pairs as well. +The higher pixel value contaminants are also potentially problematic +as they may appear to be a real source, prompting photometry to be +performed on a false object. Because of these reasons, there is a +slight preference to reject higher pixel values than lower pixel +values. + +If the number of inputs is larger than 6, then a Gaussian mixture +model analysis is run on the inputs to fit two sub populations, and +determine an the likelihood that the distribution is best described by +an uni-modal model. If this probability is less than 0.05, then the +mean is taken from the bimodal sub population with the largest +fraction of inputs, as this should exclude any sub population +comprised of high pixel value outliers. + +If this is not the case (the distribution is likely unimodal) or if +there are insufficient inputs for the mixture model analysis, the +input values are passed to an Olympic weighted mean calculation. We +set 0.2 as the fraction of the number of inputs to reject through this +process. This sets the number of bad inputs at $N_{bad} = 0.2 * +N_{input} + 0.5$, where the 0.5 term ensures at least one input is +rejected. This number is further separated into the number of low +values to exclude $N_{low} = N_{bad} / 2$, which will default to zero +if there are few inputs due to integer arithmatic, and $N_{high} = +N_{input} + N_{low} - N_{bad}$. After sorting the input values to +determine which values fall into the low and high groups, the +remaining input values have a weighted mean calculated as described +above. + +A systematic variance term is necessary to correctly scale how +discrepant points can be from the ensemble mean. If the mixture model +analysis was run, the Gaussian sigma from the largest sub population +is squared and used. If this is not available, a 0.1 scaling on the +input values is used. Each point then has a limit calculated: + +\begin{eqnarray} + limit_{mixture_model} &=& rej**2 * (variance_i + \sigma_{MM}^2) \\ + limit_{default} &=& rej**2 * (variance_i + (0.1 * value_i)**2) +\end{eqnarray} + +where $rej$ is the same factor of 4.0 used above. Each input pixel is +then compared against this limit, and the most discrepant pixel that +has $(value_i - mean)**2$ exceeding this limit is identified. If +there are suspect pixels in the set those pixels are marked for +rejection, otherwise this worst pixel is marked for rejection. +Following this, the combine and test loop is repeated for a total $0.5 +N_{input}$ iterations, or until no more pixels are rejected. + +% combineTest +%% if (Ninput > 6) { use KMM } +%% KMM: +%% Calculate KMMmu KMMsigma KMMpi KMMPunimodal +%% SumWeights = sum(pixelWeights) +%% SysVar = KMMSigma**2 OR (sys * pixelData[i])**2 +%% pixelLimts[i] = rej**2 * (pixelVariances[i] + sysVar) +% Iterate 0.5 * Ninput times (at least once) +%% Ninput = 1 => accept +%% Ninput = 2 => if (0.5 * (A - B))**2 > rej**2 * (pixelVariance[A] + pixelVariance[B] + (sys * A)**2 + (sys * B)**2) +%% then if (suspect) mark reject else mark inspect +%% Else => if (useKMM and Punimodal < 0.05) median = KMMmean +%% => else median = combinationWeightedOlympic{} +%% => if (pixelData - median)**2 > pixelLimits[i] then find single worst deviant pixel value +%% then => if suspect (mark reject) else (mark reject worst deviant pixel value) + + +%% combinationWeightedOlympic => +%% numBad = frac * Ninput + 0.5 +%% low = numBad / 2, high = low + numGood - numBad +%% sort(values) => +%% if (i > low && i <= high) { sumValues = sum_i(values * weights); sumWeight = sum_i(weights) +%% return (sumValues / sumWeight) + +% obtain lists of inspect and reject pixels. + +% normalize:? +% float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation +% psBinaryOp(ro->image, ro->image, ``*'', psScalarAlloc(normalise, PS_TYPE_F32)); +% psBinaryOp(ro->variance, ro->variance, ``*'', psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32)); + +With the initial list of rejected pixels generated, a rejection mask +is made by constructing an empty image that has the rejected pixels +set to a value of 1.0. This image is then convolved with a 5 pixel +FWHM 0-order ISIS kernel. Any pixels that are above the threshold of +0.5 are marked as bad and will be rejected in the final convolution. +If more than 10\% of all pixels from an input image are rejected, then +that entire image is rejected as well. + +% PIXEL REJECTION +% Construct 15-pixel wide ISIS kernel with 5 pixel FWHM 0-order. +% Construct image of pixels to inspect and convolve with kernel (normalize out kernel power) +% Determine pixels are bad if they're larger than THRESHOLD.MASK = 0.5. +% If more than IMAGE.REJ = 0.1 fraction of pixels are rejected, the entire image is rejected. + + +\czwdraft{I'm not entirely sure why we do what appears to be a similar + operation twice. It also seems odd that this is in the CombineFinal + step, and not in the Reject step.} Finally, the rejected pixels are +allowed to grow to include pixels that are neighbors to many rejected +pixels. The ISIS kernel used in the previous step is used to +determine the largest square box that contains under the limit of +$0.25 * \sum_{x,y} kernel**2$. This box is then convolved with the +rejected pixel mask to reject their neighbors. + +This final list of rejected pixels is passed to the final combination +pass, which does not iterate, and simply excludes the rejected +pixels. \czwdraft{This is a bad paragraph.} + +\czwdraft{We make the stacked image, the stacked variance, the stack + mask, the exposure time mask, the exp weight containing the weighted + exposure times, and a number image, containing the number of inputs + used for each pixel.} + + + + + +% FINAL COMBINE +% Grow rejected pixels +%% set threshold of (POOR.FRACTION = 0.25) * sum(kernel)**2 +%% Choose the largest square box that contains just under that threshold. +%% Convolve that box with the rejected pixels to grow them. +% Run combination pass again, but without doing rejection, simply applying the rejection lists already calculated. +% :: +% if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) { +% iter = 0 +% combineRej = NAN +% combineSys = NAN +% combineDiscard = NAN +% if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej, +% combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) { + + +The convolved stack products are not retained, as the convolution +reduces the resolution of the final image. Instead, we apply the +normalizations and rejected pixel maps generated from the convolved +stack process to the original unconvolved input images. This produces +an unconvolved stack that has the optimum image quality possible from +the input images. Not convolving does mean that the PSF shape changes +somewhat across the image, as the different FWHM of the input images +print through in the different regions in which they have contributed +to the final image. + +% UNCONVOLVED IMAGE +% if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false, true)) { +% no grow + +% only retain unconvolved products. + + +One benefit of producing the final stacked image from the weighted +mean of the unrejected input images is that faint sources do not have +their contribution removed as much as median filtering would allow. +\czwdraft{I did something to prove this once, but can't find it right + now. Comparing the ppStack output catalog to one constructed from a + simple median filtered stack shows that the ppStack catalog detects + sources up to 0.XX magnitudes fainter than the median stack. This + does increase the possibility of false positives.} @@ -755,2 +1327,5 @@ \end{document} + + +% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation