Index: trunk/doc/release.2015/ps1.detrend/detrend.tex =================================================================== --- trunk/doc/release.2015/ps1.detrend/detrend.tex (revision 39843) +++ trunk/doc/release.2015/ps1.detrend/detrend.tex (revision 39844) @@ -36,5 +36,4 @@ \newcommand{\ippstage}[1]{\textsc{#1}} \newcommand{\asinh}{\mathop{\rm asinh}\nolimits} - % Pick a terse version of the title here; @@ -127,8 +126,18 @@ reduction of the Pan-STARRS archival data. The first two reductions were used internally for pipeline optimization and the development of -the initial photometric and astrometric reference catalog. The +the initial photometric and astrometric reference catalog \citep{ps1_reference_catalog}. The products from these reductions were not publicly released, but have been used to produce a wide range of scientific papers from the Pan-STARRS 1 Science Consortium members. + +\czwdraft{Nigel: you mention calibrating to the reference catalog without telling us +what this is composed of (maybe this is in a different section, but would be +nice to have some idea here).} + +\czwdraft{Can we get around this point by simply adding a reference to + the paper describing the reference catalog? It's not really part of + the detrending process, and is discussed here mostly to give an + overview of the stages, and later to find sources of ghosts for + masking.} The Pan-STARRS image processing pipeline (IPP) is described elsewhere @@ -179,6 +188,4 @@ 24 hours of the initial set of observations \citep{WainscoatXXX}. -\czwdraft{Should there be a discussion of any header keywords/OTA file formats?} - Section \ref{sec:detrending} provides an overview of the detrending process that corrects the instrumental signatures of GPC1, with @@ -193,7 +200,5 @@ \ref{sec:discussion}. - -\czwdraft{Is this a sufficient explanation? Also, is this the right - place for it?} Image products presented in figures have been +Image products presented in figures have been mosaicked to arrange pixels as follows. Single cell images are arranged such that pixel $(1,1)$ is at the lower left corner. Images @@ -222,4 +227,9 @@ \label{sec:detrending} +\czwdraft{Nigel: I forgot: when we are talking about the various bias corrections it might be +worth pointing out that we expect these to be more of an issue in the g-band +(and maybe r?) where read noise is a significant contributor. +} + Ensuring a consistent and uniform detector response across the three-degree diameter field of view of the GPC1 camera is essential to @@ -229,5 +239,5 @@ dependent detector glows, and flat field correction to remove pixel to pixel response functions. We also construct fringe correction for the -reddest data in the y filter, to remove the interference patterns that +reddest data in the \yps{} filter, to remove the interference patterns that arise in that filter due to the variations in the thickness of the detector surface. @@ -349,5 +359,5 @@ \end{minipage} - \caption{Example of a profile cut along the y-axis through a bright star on exposure o5677g0123o OTA11 in cell xy60 (left panel) and on the subsequent exposure o5677g0124o (right panel). In both figures, the green points show the image corrected with all appropriate detrending steps, but without burntool applied, illustrating the amplitude of the persistence trails. The red points show the same data after the burntool correction, which reduces the impact of these features. Both exposures are in the g-filter with exposure times of 43s} + \caption{Example of a profile cut along the y-axis through a bright star on exposure o5677g0123o OTA11 in cell xy60 (left panel) and on the subsequent exposure o5677g0124o (right panel). In both figures, the green points show the image corrected with all appropriate detrending steps, but without burntool applied, illustrating the amplitude of the persistence trails. The red points show the same data after the burntool correction, which reduces the impact of these features. Both exposures are in the \gps{} filter with exposure times of 43s} \end{figure} @@ -564,5 +574,5 @@ % \end{subfigure} \end{minipage} - \caption{An example of the dark model application to exposure o5677g0123o, OTA23 (2011-04-26, 43s g-filter). The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, and the detector non-linearity corrected. The right panel, shows the same exposure with the dark applied in addition to the processing shown on the left.} + \caption{An example of the dark model application to exposure o5677g0123o, OTA23 (2011-04-26, 43s \gps{} filter). The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, and the detector non-linearity corrected. The right panel, shows the same exposure with the dark applied in addition to the processing shown on the left.} \end{figure} @@ -570,5 +580,5 @@ \centering \includegraphics[width=0.9\hsize,angle=0,clip]{images/B_profile_ex.png} - \caption{Example showing a profile cut across exposure o5676g0195, OTA67 (2011-04-25, 43s g-filter). The entire first row of cells (xy00-xy07) have had a median calculated along each pixel column on the OTA mosaicked image. Arbitrary offsets have been applied to shift the curves to not overlap. The top curve (in purple) shows the initial raw profile, with no dark model applied. The next curve (in green) shows the smoother profile after applying the correct B-mode dark model. Applying the incorrect A-mode dark results in the blue curve, which shows a significant increase in gradients across the cells. The orange curve shows the result of the PATTERN.CONTINUITY correction. Although this creates a larger gradient across the mosaicked images, it decreases the cell-to-cell level changes. The final yellow curve shows the final image profile after all detrending and background subtraction, and has not had an offset applied. The bright source at the cell xy00 to xy01 transition is a result of a large optical ghost, which due to the area covered, increases the median level more than the field stars.} + \caption{Example showing a profile cut across exposure o5676g0195, OTA67 (2011-04-25, 43s \gps{} filter). The entire first row of cells (xy00-xy07) have had a median calculated along each pixel column on the OTA mosaicked image. Arbitrary offsets have been applied to shift the curves to not overlap. The top curve (in purple) shows the initial raw profile, with no dark model applied. The next curve (in green) shows the smoother profile after applying the correct B-mode dark model. Applying the incorrect A-mode dark results in the blue curve, which shows a significant increase in gradients across the cells. The orange curve shows the result of the PATTERN.CONTINUITY correction. Although this creates a larger gradient across the mosaicked images, it decreases the cell-to-cell level changes. The final yellow curve shows the final image profile after all detrending and background subtraction, and has not had an offset applied. The bright source at the cell xy00 to xy01 transition is a result of a large optical ghost, which due to the area covered, increases the median level more than the field stars.} \label{fig:dark switching} \end{figure} @@ -627,5 +637,5 @@ % \end{subfigure} \end{minipage} - \caption{An example of the video dark model application to exposure o5677g0123o, OTA22 (2011-04-26, 43s g-filter), which has a video cell located in cell xy16. The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, the detector non-linearity corrected, and a regular dark applied. The right panel, shows the same exposure with a video dark applied instead of the standard dark. The main impact of this change is the improved correction of the corner glows, which are oversubtracted with the standard dark.} + \caption{An example of the video dark model application to exposure o5677g0123o, OTA22 (2011-04-26, 43s \gps{} filter), which has a video cell located in cell xy16. The left panel shows the image data mosaicked to the OTA level, and has had the static mask applied, the overscan subtracted, the detector non-linearity corrected, and a regular dark applied. The right panel, shows the same exposure with a video dark applied instead of the standard dark. The main impact of this change is the improved correction of the corner glows, which are oversubtracted with the standard dark.} \label{fig:video_darks} \end{figure} @@ -642,5 +652,7 @@ noise to increase as the row is read out. As a result of this increased noise, more sources are detected in the higher noise regions -when the read noise is assumed constant across the readout. To +when the read noise is assumed constant across the readout. Read noise is the + +To mitigate this noise gradient, we constructed an initial set of noisemap images by measuring the median variance on bias frames. The @@ -743,9 +755,18 @@ The PATTERN.ROW correction is used to remove any remaining row-by-row -bias variation, and the PATTERN.CELL and PATTERN.CONTINUITY -corrections attempt to ensure that the cells of a given OTA are -consistent with the other cells on that OTA. +bias variation, and the PATTERN.CONTINUITY correction attempts to +ensure that the cells of a given OTA are consistent with the other +cells on that OTA. \subsubsection{Pattern Row} +%% Statistics so I have them written down somewhere +%% chipProcessedImfile.bg/bg_stdev by filter for XY33 (a ``good'' chip) +%% filter bg_mean stdev median Qsig bg_stdev_mean stdev median Qsig +%% g 36.37422026669 64.64175104057 32.693 6.10284 14.696938349131 78.80460307171 8.8401 0.5417843 +%% r 200.96143304525 471.87743546238 117.105 94.55608 33.854672792146 79.01642728089 13.4564 5.3771355 +%% i 447.00504994458 938.38517801037 286.810 154.71397 57.298335510188 99.38392923935 20.0217 24.2254723 +%% z 317.54933679054 390.38930252748 241.014 114.13316 48.359069000176 94.44452756094 17.9404 9.1535209 +%% y 371.09019536218 293.57439970375 288.481 133.38769 43.724342221691 135.04286534327 19.9029 7.5396461 + % http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study As discussed above in the dark and noisemap sections, certain @@ -772,4 +793,9 @@ the sky. +These row-by-row variations have the largest impact on data taken in +the \gps{} filter, as the read noise is the dominant noise source in +that filter. At longer wavelengths, the noise from the Poissonian +variation in the sky level increases. Although the PATTERN.ROW correction is still applied to data taken in the other filters, + This correction was required on all cells on all OTAs prior to 2009-12-01, at which point a modification of the camera electronics @@ -841,10 +867,11 @@ % \end{subfigure} \end{minipage} - \caption{Example of the PATTERN.ROW correction on exposure o5379g0103o OTA57 cell xy00 (i-filter 45s). The left panel shows the cell with all appropriate detrending except the PATTERN.ROW, and the right shows the same cell with PATTERN.ROW applied. The correction reduces the correlated noise on the right side, which is most distant from the read out amplifier. There is a slight over subtraction along the rows near the bright star.} + \caption{Example of the PATTERN.ROW correction on exposure o5379g0103o OTA57 cell xy00 (\ips{} filter 45s). The left panel shows the cell with all appropriate detrending except the PATTERN.ROW, and the right shows the same cell with PATTERN.ROW applied. The correction reduces the correlated noise on the right side, which is most distant from the read out amplifier. There is a slight over subtraction along the rows near the bright star.} \end{figure} \subsubsection{Pattern Continuity} -As the PATTERN.CELL correction was insufficient in many situations, we +After previous attempts to ensure that adjacent cells on an OTA +matched background levels were insufficient in many situations, we designed a replacement correction that would reduce the background distortion for large objects. In addition, studies of the background @@ -855,19 +882,16 @@ horizontally across an OTA, and as the background model fits a smooth sky level, this induces over and under subtraction at the cell -boundaries. As the PATTERN.CELL was designed to correct changes only -in the median value between cells, it could not adequately resolve -this higher order issue. - -The replacement for PATTERN.CELL is the PATTERN.CONTINUITY correction, -which attempts to match the edges of a cell to those of its neighbors. -For each cell, a thin box 10 pixels wide on each edge is extracted and -the median value of unmasked values calculated for that box. These -median values are then used to construct a vector of differences -$\Delta_i = \sum_{j} Edge_{i} - Edge_{j}$, along with a matrix of -associations $A_{i,i'} = \sum_{j} \delta(i,j) \delta(j,i')$ denoting -which cell boundaries are adjacent. By solving the system $A x = -diff$, we find the set of offsets $x_i$ to be applied to each cell to -ensure the minimum differences between all cell edges and their -neighbors. +boundaries. + +The PATTERN.CONTINUITY correction, attempts to match the edges of a +cell to those of its neighbors. For each cell, a thin box 10 pixels +wide on each edge is extracted and the median value of unmasked values +calculated for that box. These median values are then used to +construct a vector of differences $\Delta_i = \sum_{j} \mathrm{Edge}_{i} - +\mathrm{Edge}_{j}$, along with a matrix of associations $A_{i,i'} = \sum_{j} +\delta(i,j) \delta(j,i')$ denoting which cell boundaries are adjacent. +By solving the system $A x = \Delta$, we find the set of offsets $x_i$ +to be applied to each cell to ensure the minimum differences between +all cell edges and their neighbors. For OTAs that initially show the saw tooth pattern, the effect of this @@ -895,7 +919,7 @@ wavelength of the light becomes comparable to the thickness of the detectors. Visually inspecting the images shows that the fringing is -most prevalent in the y filter images, with negligible fringing in the +most prevalent in the \yps{} filter images, with negligible fringing in the other bands. As a result of this, we only apply a fringe correction -to the y filter data. +to the \yps{} filter data. The fringe used for PV3 processing was constructed from a set of 20 @@ -925,15 +949,10 @@ \centering \begin{minipage}{0.5\hsize} - \includegraphics[width=1.0\hsize,angle=0,clip]{images/o5220g0025o_XY53_nofringe.png} -% \caption{(a)} -% \end{subfigure}% -% \begin{subfigure}[]{.45\hsize} + \includegraphics[width=1.5\hsize,angle=0,clip]{images/o5220g0025o_XY53_nofringe.png} \end{minipage}% \begin{minipage}{0.5\hsize} - \includegraphics[width=1.0\hsize,angle=0,clip]{images/o5220g0025o_XY53_fringe.png} -% \caption{(b)} -% \end{subfigure} + \includegraphics[width=1.5\hsize,angle=0,clip]{images/o5220g0025o_XY53_fringe.png} \end{minipage} - \caption{Example of the y-filter fringe pattern on exposure o5220g0025o OTA53 (y-filter 30s). The left panel shows the OTA mosaic with all detrending except the fringe correction, while the right shows the same including the fringe correction. Both images have been smoothed with a Gaussian with $\sigma = 3$ pixels to highlight the faint and large scale fringe patterns. \czwdraft{See if there's a way to have mana produce images larger than the screen size.}} + \caption{Example of the \yps{} filter fringe pattern on exposure o5220g0025o OTA53 (\yps{} filter 30s). The left panel shows the OTA mosaic with all detrending except the fringe correction, while the right shows the same including the fringe correction. Both images have been smoothed with a Gaussian with $\sigma = 3$ pixels to highlight the faint and large scale fringe patterns. \czwdraft{See if there's a way to have mana produce images larger than the screen size.}} \label{fig: fringe example} \end{figure} @@ -983,5 +1002,5 @@ The final step of mask construction is to examine the detector for bright columns and other static pixel issues. This is first done by -processing a set of 100 i filter science images in the same fashion as +processing a set of 100 \ips{} filter science images in the same fashion as for the DARKMASK. A median image is constructed from these inputs along with the per-pixel variance. These images are used to identify @@ -1124,5 +1143,5 @@ Due to imperfections in the anti-reflective coating on the optical surfaces of GPC1, bright sources can also result in large out of focus -objects, particularly in the g-filter data. These objects are the +objects, particularly in the \gps{} filter data. These objects are the result of light reflecting back off the surface of the detector, reflecting again off the lower surfaces of the optics (particularly @@ -1182,10 +1201,10 @@ \tablehead{\colhead{Filter}&\colhead{$m_{inst}$}} \startdata - g & -16.5 \\ - r & -20.0 \\ - i & -25.0 \\ - z & -25.0 \\ - y & -25.0 \\ - w & -20.0 \\ + \gps{} & -16.5 \\ + \rps{} & -20.0 \\ + \ips{} & -25.0 \\ + \zps{} & -25.0 \\ + \yps{} & -25.0 \\ + \wps{} & -20.0 \\ \enddata \label{tab:ghost_magnitudes} @@ -1196,5 +1215,5 @@ \centering \includegraphics[width=0.9\hsize,angle=0,clip]{images/full_fpa_ghosts.jpg} - \caption{Example of the full GPC1 field of view illustrating the sources and destinations of optical ghosts on exposure o5677g0123o (2011-04-26, 43s g-filter). The bright stars on OTA33 and OTA44 result in nearly circular ghosts on the opposite OTA. In contrast, the trio of stars on OTA11 result in very elongated ghosts on OTA66.} + \caption{Example of the full GPC1 field of view illustrating the sources and destinations of optical ghosts on exposure o5677g0123o (2011-04-26, 43s \gps{} filter). The bright stars on OTA33 and OTA44 result in nearly circular ghosts on the opposite OTA. In contrast, the trio of stars on OTA11 result in very elongated ghosts on OTA66.} \end{figure} @@ -1202,11 +1221,41 @@ \label{sec:glints} -Prior to \czwdraft{DATE}, a reflective surface at the edge of the -camera aperture was incompletely screened to light passing through the +% I finally tracked it down: +%% > On 8/26/2010 9:24 AM, John Tonry wrote: +%% > +%% > Gene, +%% > +%% > This is a bit of a case of the dog that didn't bark, but the shutter mask +%% > went in on Tuesday. +%% > +%% > Can you can tell us whether +%% > +%% > a) it's helped the shutter glint problem and +%% > b) whether there's any discernable vignetting anywhere? +%% > +%% > - John + +%% On Thu, Aug 26, 2010 at 4:00 PM, Chris Waters wrote: + +%% > I'm not entirely sure why I'm not on the ps-ipp mailing list, but +%% > Heather forwarded this to me. I compared 240 exposures from +%% > 2010-08-22/ThreePi/y.00000 and 2010-08-25/ThreePi/y.00000. +%% > +%% > a) For the 22nd, I counted 28 star glints visible. For the 25th, I +%% > counted maybe 0-2 (I think the first is a conveniently placed satellite, +%% > and the other has a companion, so I think it's actually a moon glint). +%% > +%% > b) I was going to compare flat field images, but we don't have any +%% > from after the mask was applied. Blinking between a few pairs of the +%% > 240x2 exposures does not show any vignetting that I can detect from +%% > the IPP jpeg mosaics. + +Prior to 2010-08-24, a reflective surface at the edge of the camera +aperture was incompletely screened to light passing through the telescope. Sources brighter than $m_{inst} = -21$ that fell on this reflective surface resulted in light being scattered across the detector surface in a long narrow glint. This surface was physically -masked on \czwdraft{DATE}, removing the possibility of glints in -subsequent data, but that taken prior have a dynamic mask constructed +masked on 2010-08-24, removing the possibility of glints in subsequent +data, but that taken prior have an advisory dynamic mask constructed when a reference source falls on the focal plane within one degree of the detector edge. This mask is 150 pixels wide, with length $L = @@ -1244,5 +1293,5 @@ \centering \includegraphics[width=0.9\hsize,angle=0,clip]{images/glint_example_o5379g0103o.jpg} - \caption{Example of a glint on exposure o5379g0103o (2010-07-02, 45s i-filter). The source star out of the field of view creates a long reflection that extends through OTA73 and OTA63.} + \caption{Example of a glint on exposure o5379g0103o (2010-07-02, 45s \ips{} filter). The source star out of the field of view creates a long reflection that extends through OTA73 and OTA63.} \end{figure} @@ -1271,5 +1320,5 @@ \centering \includegraphics[width=0.9\hsize,angle=0,clip]{images/o6802g0338o_XY51_b1.jpg} - \caption{Example of saturated star, with diffraction spikes extending from the core on exposure o6802g0338o, OTA51 (2014-05-25, 45s g-filter).} + \caption{Example of saturated star, with diffraction spikes extending from the core on exposure o6802g0338o, OTA51 (2014-05-25, 45s \gps{} filter).} \label{fig:saturated star} \end{figure} @@ -1295,5 +1344,11 @@ calculations to estimate the masking fraction. 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. \czwdraft{I need a percentage here.} +arcseconds per pixel, translates to a 20930 FPA pixel radius. Summing +mask fractions from these three contributions within the unvignetted +field of view results in an average of $\sim 20\%$ masking fraction +across the field of view. Dynamic masking adds an additional $2-3\%$ +on average, with advisory burntool masking contributing the largest +single component. + %% 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; @@ -1313,27 +1368,131 @@ %% | 0.21130344126869 | 0.00013634812877977 | 0.02163070300815 | -Summing mask fractions from these three contributions within the -unvignetted field of view results in an average of $\sim 20\%$ masking -fraction across the field of view. Dynamic masking adds an additional -$2-3\%$ on average, with advisory burntool masking contributing the -largest single component. \subsection{Background subtraction} \label{sec:background} + +\czwdraft{Nigel: 2.10 The background section is rather short, given all the fuss DRAVG made +about it. What is done to eliminate contamination by bright objects - isn't +there some sort of clipping? We also have a confusing number of ``bins'' in the +text (``These bins have 10000 .... a binned cumulative distribution is +generated. These bins are interpolated ... levels. Repeating this across all +bins ...''). There is no mention of the fact that this will subtract real +astrophysics backgrounds if they are on a suitably large scale, or of the fact +that the subtraction is not perfect (don't I remember that the stacks end up +with a non-zero background which scales with the number of input warps?). +} + +\czwdraft{Based on the wiki page on 2014-05-21 the stack background issue should be resolved.} Once all other detrending is done, the pixels from each cell are mosaicked into the full $4846\times{}4868$ pixel OTA image. A background model for the full OTA is then determined prior to the -photometric analysis. The mosaicked image is binned into -$800\times{}800$ pixel bins, centered on the image center, and -overlapping by a factor of 2 in both axes. These bins have 10000 -random samples drawn, and a binned cumulative distribution function is -generated. These bins are interpolated to find the best mean value at -the $50\%$ level, as well as the distribution $\sigma$ by estimating -from the $32\%$ and $68\%$ levels. Repeating this across all bins -results in a $13\times{}13$ grid of background bins, which are -bilinearly interpolated to generate the background model to subtract. -Each object in the photometric catalog has a SKY and SKY\_SIGMA value -based on this model as well. +photometric analysis. The mosaicked image is subdivided into +$800\times{}800$ pixel segments that define each pixel of the +background model, with the segments centered on the image center, and +overlapping adjacent subdivisions by 400 pixels. These overlaps help +smooth the background model, as adjacent model pixels share input +pixels. + +From each subdivision, 10000 random unmasked pixels are drawn. In the +case where the mask fraction is large (such as on OTAs near the edge +of the field of view), and there are insufficient unmasked pixels to +meet this criterion, all possible unmasked pixels are used instead. +If this number is still small (less than 100 good pixels), the +subdivision does not have a background model calculated, and instead, +the value assigned to that model pixel is set as the average of the +adjacent model pixels. This allows up to eight neighboring background +values to be used to patch these bad pixels. + +For the remaining subdivisions that have sufficient unmasked pixels +for the background to be measured, the pixel values are used to +calculate a set of robust statistics for the initial background guess. +The minimum and maximum of the values are found, and checked to ensure +that these are not the same value, which would indicate some problem +with the input values. The values are then inserted into a histogram +with 1000 bins between the minimum and maximum values, and again +checked for issues with the inputs by ensuring that the bin with the +most input pixels does not contain more than half of the input values. +In this case, the minimum and maximum do not constrain the true +distribution of the input values well, and any values outside of the +20 bins closest to the bin with the peak are masked for future +consideration. A cumulative distribution is then constructed from the +histogram, which saves the computational cost of sorting all the input +values. The bins containing the 50-percentile point, as well as the +15.8\%, 84.1\% ($\pm 1 \sigma$), 30.8\%, 69.1\% ($\pm 0.5 \sigma$), +2.2\%, and 97.7\% ($\pm 2 \sigma$) points are identified in this +cumulative histogram. These bins, and the two bins to either side are +then linearly interpolated to identify the pixel value corresponding +to these points in the distribution. The 50\% point is set as the +median of the pixel distribution, with the standard deviation of the +distribution set as the median of the $\sigma$ values calculated from +the $0.5 * (\sigma_{+1} - \sigma_{-1})$, $\sigma_{+0.5} - +\sigma_{-0.5}$, and $0.25 * (\sigma_{+2} - \sigma_{-2})$ differences. +If this measured standard deviation is smaller than 3 times the bin +size, then all points more than 25 bins away from the calculated +median are masked, and the process is repeated until the bin size is +sufficiently small to ensure that the distribution width is well +sampled. Once this iterative process converges, or 20 iterations are +run, the 25- and 75-percentile values are found by interpolating the 5 +bins around the expected bin as well, and the count of the number of +input values within this inner 50-percentile region, $N_{50}$ is +calculated. + +These initial statistics are then used as the starting guesses for a +second calculation of the background level that attempts to fit the +distribution with a Gaussian. All pixels that were masked in the +initial calculation are unmasked, and a histogram is again constructed +of the values, with a binsize set to $\sigma_{guess} / \left( N_{50} / +500 \right)$. With this binsize, we expect that a bin at $\pm 2 +\sigma$ will have approximately 50 input points, which gives a +Poissonian signal to noise estimate around 7. In the case where +$N_{50}$ is small (due to a poorly populated input image), this bin +size is fixed to be no larger than the guess of the standard +deviation. The endpoints of the histogram are clipped based on the +input guesses, such that any input point with a value more than $5 +\sigma_{guess}$ away from the input mean are excluded from +consideration. + +Two second order polynomial fits are then performed to the logarithm +of the histogram counts set at the midpoint of each bin. The first +fit considers the ``lower half'' of the distribution, under the +assumption that deviations from a normal distribution are caused by +real astrophysical sources that will be brighter than the true +background level. From the bin with most pixel values, the lower +bound is set by searching for the first bin from the peak that has +fewer inputs than 25\% of the peak. A similar search is performed for +the upper bound, but with a criterion that the bin has fewer than 50\% +of the peak. On both sides of the peak, the bounds are adjusted to +ensure that at least seven bins, equally distributed around the peak, +are used. The second fit is symmetric, fitting both sides of the +distribution out to the point where the bin contains fewer than 15\% +of the peak value. The same seven-bin constraint is used for this +fit. The Gaussian mean and standard deviation are calculated from the +polynomial coefficients, and the symmetric fit results are accepted +unless the lower-half fit results in a smaller mean. This process is +repeated again if the calculated standard deviation is not larger than +75\% of the initial guess (suggesting an issue with the initial bin +size). + +With this two-stage calculation performed across all subdivisions of +the mosaicked OTA image, and missing model pixels filled with the +average of their neighbors, the final background model is stored on +disk as a $13\times{}13$ image with header entries listing the binning +used. The full scale background image is then constructed by +binlinearly interpolating this binned model, and this is subtracted +from the science image. Each object in the photometric catalog has a +SKY and SKY\_SIGMA value that is the evaluation of this model at the +location of that object. + +Although this background modeling process works well for most of the +sky, astronomical sources that are large compared to the +$800\times{}800$ pixel subdivisions can bias the calculated background +level high, resulting in an oversubtraction near that object. The +most common source that can cause this issue are large galaxies, which +can have their own features modeled as being part of the background. +For the specialized processing of M31, which covers an entire pointing +of GPC1, the measured background was added back to the \ippstage{chip} +stage images, but this special processing was not used for the large +scale $3\Pi$ PV3 reduction. %% * Magic @@ -1452,9 +1611,10 @@ & 964 & 2010-09-01 00:00:00 & 2011-05-01 00:00:00 & \\ & 965 & 2011-05-01 00:00:00 & & \\ - FLAT & 300 & 2009-12-09 00:00:00 & & g filter \\ - & 301 & 2009-12-09 00:00:00 & & r filter \\ - & 302 & 2009-12-09 00:00:00 & & i filter \\ - & 303 & 2009-12-09 00:00:00 & & z filter \\ - & 304 & 2009-12-09 00:00:00 & & y filter \\ + FLAT & 300 & 2009-12-09 00:00:00 & & \gps{} filter \\ + & 301 & 2009-12-09 00:00:00 & & \rps{} filter \\ + & 302 & 2009-12-09 00:00:00 & & \ips{} filter \\ + & 303 & 2009-12-09 00:00:00 & & \zps{} filter \\ + & 304 & 2009-12-09 00:00:00 & & \yps{} filter \\ + & 305 & 2009-12-09 00:00:00 & & \wps{} filter \\ FRINGE & 296 & 2009-12-09 00:00:00 & & \\ ASTROM & 1064 & 2008-05-06 00:00:00 & & \\ @@ -1530,5 +1690,56 @@ 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?} +parameters for the locally linear grid. These parameters are stored +in a string listing the reference position in the chip coordinate +frame, the slope of the relation in the warp x axis, and the slope of +the relation in the warp y axis. From these keywords, any position in +the warp can be mapped back to the location in any of the input OTA +images. + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_sci.jpg} + \caption{Example of the warp image for skycell skycell.2047.005 + centered at ($\alpha,\delta$) = (179.763, 32.1899) for exposure + o4985g0073o, (2009-06-03, 30s \zps{} filter). The data from six + OTAs contribute to this image, although they are all truncated by + the skycell boundaries. This skycell image is aligned such that + north points to the top of the image, and east to the left. The + contributing OTAs are from the right half of the detector, with + OTA24 contributing the most pixels, and originally have the + positive y axis pointing to the southwest in this warped image, + with the positive x axis to the northwest.} + \label{fig:warp image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_wt.jpg} + \caption{Example of the warp variance image for skycell + skycell.2047.005 of exposure o4985g0073o, the same as in Figure + \ref{fig:warp image}. This variance map retains information about + the higher flux levels that were found in burntool corrected + persistence trails, which appear here as streaks along the + original OTA y axis. The amplifier glows that are corrected in + the dark model are also more visible in the corners of the cells + in OTA24. As both of these effects are corrected in the science + image, there are no significant features visible there.} + \label{fig:warp variance} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/warp_1046511_sci.jpg} + \caption{Example of the warp mask image for skycell skycell.2047.005 + of exposure o4985g0073o, the same as in Figure \ref{fig:warp + image}. This mask image shows the many small defects removed + from the image, along with larger advisory trails on corrected + burntool trails. The saturated cores of the bright stars are also + masked, along with the diffraction spikes found on these stars. + In addition OTA24 shows the precautionary crosstalk bleed masks + for the two brightest stars applied to all cells within the same + row.} +\end{figure} + % Read all images and astrometry @@ -1558,10 +1769,10 @@ system, they can then be combined pixel-by-pixel regardless of their original orientation. Creating a stacked image by coadding the -individual warps increases the signal to noise, allowing objects -fainter than the single image signal to noise threshold. Creating -this stack also allows a complete image to be constructed that does -not have regions masked due to the gaps between cells and OTAs. This -fully populated static sky image can also be used as a template for -subtraction to find transient sources. +individual warps increases the signal to noise, allowing for the +detection of objects that would not be sufficiently significant to be measured from a single image. +Creating this stack also allows a complete image to be +constructed that does not have regions masked due to the gaps between +cells and OTAs. This fully populated static sky image can also be +used as a template for subtraction to find transient sources. The stacked image is comprised of all warp frames for a given skycell @@ -1572,31 +1783,42 @@ 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 seeing far worse -than average, and would degrade the final output stack. For the PV3 -survey, this size represents a PSF larger than $97$th percentile in -all filters. A target PSF for the stack is constructed by finding the -maximum envelope of all input PSFs, which sets the target PSF to the -largest value among the input PSFs for a given position from the peak. -This PSF is then circularized to ensure azimuthal symmetry, which -prevents any of the input images from being deconvolved when matched -to the target. - -The input images also need to have their flux normalized to prevent +FWHM greater than 10 pixels (2.5 arcseconds), as those images have the +seeing far worse than average, and would degrade the final output +stack. For the PV3 $3\Pi$ survey, this size represents a PSF larger +than the $97$th percentile in all filters. A target PSF for the stack +is constructed by finding the maximum envelope of all input PSFs, +which sets the target PSF to the largest value among the input PSFs +for a given position from the peak. This PSF is then circularized to +ensure azimuthal symmetry, which prevents deconvolution of any of the +input images when matched to the target. + +The input images also need to have their fluxes 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 to the sum of the input exposure -times. We can determine the relative transparency for each input -image by comparing the magnitudes of matched sources between the -different images. Each image then has a normalization factor defined, -equal to $norm_{i} = (ZP_{i} - ZP_{target}) - transparency_{i} - 2.5 * -\log_{10} (t_{target} / t_{i}) - airmassTerm * (airmass_{i} - -airmass_{target})$. \czwdraft{ZP.AIRMASS is zero for all filters. - Does this simply mean that we assume any airmass differences are - folded into the transparency differences? This would simplify this - discussion quite a bit if that's the case, as we can just say that - and skip all the extra airmass terms.} +during pixel rejection. From the reference catalog 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 to the sum of +the input exposure times, regardless of if those inputs are rejected +later in the combination process. We can determine the relative +transparency for each input image by comparing the magnitudes of +matched sources between the different images. Each image then has a +normalization factor defined, equal to $\mathrm{norm}_{input} = (ZP_\mathrm{input} +- ZP_\mathrm{target}) - \mathrm{transparency}_\mathrm{input} - 2.5 * +\log_{10} (t_\mathrm{target} / t_\mathrm{input}) - +\mathrm{F}_\mathrm{airmass} * (\mathrm{airmass}_\mathrm{input} - +\mathrm{airmass}_\mathrm{target})$. For the PV3 processing, the +airmass factor $\mathrm{F}_\mathrm{airmass}$ was set to zero, such +that all flux differences from differing exposure airmasses are +assumed to be included in the zeropoint and transparency values. + + +\czwdraft{Nigel: 5. ``The ouput exposure time is set to the sum of the input exposure times.'' +True, but we should note that as warps can be rejected later on in the +stacking process this output time is notional in some sense. +Calibration - for PV3 what photometric calibration has been used at this stage +for the input warps? Should we make it clear here that pixels are not subject +to the final (any?) ubercal? +} % PREPARE @@ -1636,9 +1858,21 @@ convolution kernels can be calculated for each image. ISIS kernels \citep{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. - +pixels and polynomial orders of 6, 4, and 2. Regions around the +sources identified in the input images are extracted, convolved with +the kernel, and the residual with the target PSF used to update the +parameters of the kernel via least squares optimization. Stamps that +significantly deviate are rejected, but as the squared residual +difference will increase with increasing source flux. To mitigate +this effect, a parabola is fit to the distribution of squared +residuals as a function of source flux. Stamps that deviate from this +fit by more than $2.5\sigma$ are rejected, and not used on further +kernel fit iterations. This process is repeated twice, and the final +convolution kernel is returned. + +This convolution may change the image flux scaling, so a normalization +factor is used to correct this. This normalization factor is equal to +the ratio of $10^{-0.4 \mathrm{norm}_{input}}$ to the sum of the +kernel. The image is multiplied by this factor, and the variance by +the square of it, scaling all inputs to the common zeropoint. % MATCH @@ -1651,11 +1885,11 @@ Once the convolution kernels are defined for each image, they are used to convolve the image to match the target PSF. Any input image that -has a $\chi^2$ value greater than 4.0$\sigma$ larger than the median -value is rejected from the stack. Each image also has a weight -assigned, based on the image variance after convolution. For a given -image, the weight is equal to $W^{-1} = \langle Variance(x,y) \rangle -* f_{covariance}$, where the angle brackets denote a robust median of -the variance image, and the covariance factor $f_{covariance}$ is the -peak value of the covariance matrix of the convolution. +has a kernel match $\chi^2$ value greater than 4.0$\sigma$ larger than +the median value is rejected from the stack. Each image also has a +weight assigned, based on the image variance after convolution. A +full image weight is then calculated for each input, with the weight, +$W_\mathrm{input}$ is equal to the inverse of the median of the image +variance multiplied by the peak of the image covariance (due to the +warping process). % CONVOLVE @@ -1685,6 +1919,6 @@ \begin{eqnarray} - S_{value} &=& \sum_i\left(value_{i} * W_i\right) / \sum_i\left(W_i\right) \\ - S_{exp weight} &=& \sum_i \left(exptime_i * W_i\right) / \sum_i\left(W_i\right) \\ + \mathrm{Stack}_\mathrm{value} &=& \sum_i\left(\mathrm{value}_\mathrm{input} * W_\mathrm{input}\right) / \sum_\mathrm{inputs}\left(W_\mathrm{input}\right) \\ + \mathrm{Stack}_\mathrm{exp weight} &=& \sum_i \left(\mathrm{exptime}_\mathrm{input} * W_\mathrm{input}\right) / \sum_\mathrm{inputs}\left(W_\mathrm{input}\right) \\ \end{eqnarray} @@ -1692,5 +1926,5 @@ \begin{eqnarray} - S_{variance} &=& 1 / \sum_i \left( 1 / variance_i \right) + \mathrm{Stack}_\mathrm{variance} &=& 1 / \sum_i \left( 1 / \sigma^2_\mathrm{input} \right) \end{eqnarray} @@ -1766,12 +2000,14 @@ to reject higher pixel values than lower pixel values. -Following this initial combination, a ``testing'' loop iterates in an +Following the 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 sigmas 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. +$B$, then a check is made to see if $(0.5 * (\mathrm{value}_A - +\mathrm{value}_B))^2 > 16 * (\sigma^2_A + \sigma^2_B ++ (0.1 * \mathrm{value}_A)^2 + (0.1 * \mathrm{value}_B)^2)$, such that +the deviation of the inputs from their mean position is greater than +four times the sum of their measured uncertainties and a 10\% +systematic error term. If this is the case, neither input is trusted, +and both are flagged for rejection If the number of inputs is larger than 6, then a Gaussian mixture @@ -1787,12 +2023,13 @@ input values are passed to an Olympic weighted mean calculation. We reject $20\%$ of the number of inputs through this process. The -number of bad inputs is set to $N_{bad} = 0.2 * N_{input} + 0.5$, with -the 0.5 term ensuring 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, 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 are used in a weighted mean using the image -weights above. +number of bad inputs is set to $N_\mathrm{bad} = 0.2 * +N_\mathrm{input} + 0.5$, with the 0.5 term ensuring at least one input +is rejected. This number is further separated into the number of low +values to exclude $N_\mathrm{low} = N_\mathrm{bad} / 2$, which will +default to zero if there are few inputs, and $N_\mathrm{high} = +N_\mathrm{input} + N_\mathrm{low} - N_\mathrm{bad}$. After sorting +the input values to determine which values fall into the low and high +groups, the remaining input values are used in a weighted mean using +the image weights above. A systematic variance term is necessary to correctly scale how @@ -1804,15 +2041,16 @@ \begin{eqnarray} - limit_{mixture model} &=& 4^2 * (variance_i + \sigma_{MM}^2) \\ - limit_{default} &=& 4^2 * (variance_i + (0.1 * value_i)^2) + \mathrm{limit}_\mathrm{mixture model} &=& 4^2 * (\sigma^2_\mathrm{input} + \sigma_\mathrm{mixture model}^2) \\ + \mathrm{limit}_\mathrm{default} &=& 4^2 * (\sigma^2_\mathrm{input} + (0.1 * \mathrm{value}_\mathrm{input})^2) \end{eqnarray} 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 -until no more pixels are rejected, up to a maximum number of -iterations equal to $50\%$ of the number of inputs. +discrepant pixel that has $(\mathrm{value}_\mathrm{input} - +\mathrm{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 until no more pixels are +rejected, up to a maximum number of iterations equal to $50\%$ of the +number of inputs. % combineTest @@ -1848,10 +2086,11 @@ 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 zeroth-order ISIS kernel. Any pixels that are above the threshold of -0.5 after this mask convolution are marked as bad and will be rejected in the final combination. -If more than 10\% of all pixels from an input image are rejected, then -that entire image is rejected as well. +is made for the input warp by constructing an empty image that has the +rejected pixels from that input set to a value of 1.0. This image is +then convolved with a 5 pixel FWHM zeroth-order ISIS kernel. Any +pixels that are above the threshold of 0.5 after this mask convolution +are marked as bad and will be rejected in the final combination. If +more than 10\% of all pixels from an input image are rejected, then +the entire image is rejected as it likely has some systematic issue. % PIXEL REJECTION @@ -1862,18 +2101,16 @@ -\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 +Finally, a second pass at rejecting pixelsis conducted, by growing the +current list to include pixels that are neighbors to many rejected +pixels. The ISIS kernel used in the previous step is again 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 the neighboring pixels. This final list of -rejected pixels is passed to the final combination, which creates the -final stack values from the weighted mean of the non-rejected pixels. -Six total images are constructed for this final stack: the image, its -variance, a mask, a map of the exposure time per pixel, that exposure -time map weighted by the input image weight, and a map of the number -of inputs per pixel. +$0.25 * \sum_{x,y} kernel^2$. This square box is then convolved with +the rejected pixel mask to reject the neighboring pixels. This final +list of rejected pixels is passed to the final combination, which +creates the final stack values from the weighted mean of the +non-rejected pixels. Six total images are constructed for this final +stack: the image, its variance, a mask, a map of the exposure time per +pixel, that exposure time map weighted by the input image weight, and +a map of the number of inputs per pixel. % FINAL COMBINE @@ -1945,4 +2182,96 @@ such that: $L = \mathrm{BOFFSET} + \mathrm{BSOFTEN} \cdot \left(\exp(C / \alpha) - \exp(-C / \alpha)\right)$. + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_sci.jpg} + \caption{Example of the stack image for skycell skycell.2047.005 + centered at ($\alpha,\delta$) = (179.763, 32.1899) in the \zps{} + filter, stack\_id 3775944. This stack includes 25 input images, + including o4985g0073o the warp image in Figure \ref{fig:warp + image}, and has a combined exposure time of 870s. Combining + such a large number of input images removes the inter-cell and + inter-chip gaps, providing a fully populated image. In addition, + the combined signal allows many more faint objects to be found + than were visible on the single frame warp image.} + + \label{fig:stack image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_mask.jpg} + \caption{Example of the stack mask image for skycell + skycell.2047.005 centered at ($\alpha,\delta$) = (179.763, + 32.1899) in the \zps{} filter, stack\_id 3775944. The entire + frame is largely unmasked after combining inputs, with the only + remaining masks falling on the cores of bright stars, and in small + regions around the brighest objects where the overlapping of + diffraction spike masks have removed all inputs.} + + \label{fig:stack mask image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_wt.jpg} + \caption{Example of the stack variance image for skycell + skycell.2047.005 centered at ($\alpha,\delta$) = (179.763, + 32.1899) in the \zps{} filter, stack\_id 3775944. The variance + map for this stack is reasonably smooth, with the mottled pattern + from the inter-chip and inter-cell gaps printing through. Some + regions with higher variance are found where the number of inputs + is lower.} + + \label{fig:stack wt image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_num.jpg} + \caption{Example of the stack number image for skycell + skycell.2047.005 centered at ($\alpha,\delta$) = (179.763, + 32.1899) in the \zps{} filter, stack\_id 3775944. This map shows + the number of inputs contributing to each pixel of the output + stack. Again, the pattern of the inter-chip and inter-cell gaps + is visible, along with the mask pattern of regions with CTE + problems (visible in the upper right corner). } + + \label{fig:stack num image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_exp.jpg} + \caption{Example of the stack exposure time image for skycell + skycell.2047.005 centered at ($\alpha,\delta$) = (179.763, + 32.1899) in the \zps{} filter, stack\_id 3775944. As all input + warps had the same 30s exposure time, this map essentially + recreates the number map, with units of seconds of exposure + instead of number of inputs contributing to a given pixel.} + + \label{fig:stack exp image} +\end{figure} + +\begin{figure} + \centering + \includegraphics[width=0.9\hsize,angle=0,clip]{images/stack_3775944_expwt.jpg} + \caption{Example of the stack weighted exposure image for skycell + skycell.2047.005 centered at ($\alpha,\delta$) = (179.763, + 32.1899) in the \zps{} filter, stack\_id 3775944. This map shows + the weighted average exposure time, as described in the text. It + is similar to the simple exposure time map, but shows how some + input exposures have their contributions weighted down due to the + observed larger image variances.} + + + \label{fig:stack exp wtimage} +\end{figure} + + + + + + \section{Discussion}