Superresolution Reconstruction of Severely Undersampled Point-spread Functions Using Point-source Stacking and Deconvolution
Teresa Symons, Michael Zemcov, James Bock, Yun-Ting Cheng, Brendan Crill, Christopher Hirata, Stephanie Venuto
DDraft version February 3, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Superresolution Reconstruction of Severely Undersampled Point-spread Functions Using Point-sourceStacking and Deconvolution
Teresa Symons , Michael Zemcov ,
1, 2
James Bock,
2, 3
Yun-Ting Cheng , Brendan Crill , Christopher Hirata, and Stephanie Venuto Center for Detectors, Rochester Institute of Technology, 1 Lomb Memorial Drive, Rochester, NY 14623, USA Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA (Accepted for publication in ApJS)
ABSTRACTPoint-spread function (PSF) estimation in spatially undersampled images is challenging becauselarge pixels average fine-scale spatial information. This is problematic when fine-resolution details arenecessary, as in optimal photometry where knowledge of the illumination pattern beyond the nativespatial resolution of the image may be required. Here, we introduce a method of PSF reconstructionwhere point sources are artificially sampled beyond the native resolution of an image and combinedtogether via stacking to return a finely sampled estimate of the PSF. This estimate is then deconvolvedfrom the pixel-gridding function to return a superresolution kernel that can be used for optimallyweighted photometry. We benchmark against the <
1% photometric error requirement of the upcomingSPHEREx mission to assess performance in a concrete example. We find that standard methods likeRichardson–Lucy deconvolution are not sufficient to achieve this stringent requirement. We investigatea more advanced method with significant heritage in image analysis called iterative back-projection(IBP) and demonstrate it using idealized Gaussian cases and simulated SPHEREx images. In testingthis method on real images recorded by the LORRI instrument on New Horizons, we are able to identifysystematic pointing drift. Our IBP-derived PSF kernels allow photometric accuracy significantly betterthan the requirement in individual SPHEREx exposures. This PSF reconstruction method is broadlyapplicable to a variety of problems and combines computationally simple techniques in a way that isrobust to complicating factors such as severe undersampling, spatially complex PSFs, noise, crowdedfields, or limited source numbers. INTRODUCTIONWhether explicitly or implicitly, astrophysicists re-quire knowledge of how the brightness in an image re-lates to true brightness on the sky to interpret the shapesand intensities of observed sources of emission. The in-strument response function of a telescope is the transferfunction that maps between quantities measured at thefocal plane to the physical intensity on the sky. Fornatively imaging array detectors, or images constructedfrom time-domain scans of the sky, this information isoften called the point-spread function (PSF), which isdefined to be the response of a focused imaging sys-
Corresponding author: Teresa [email protected] tem to a point source (see Gai & Cancelliere 2007 fora review). PSF estimation is a foundational problem inastronomical image analysis and interpretation.Measuring the PSF of an instrument can be challeng-ing. For some instruments, it is possible to measure thePSF using collimated images of unresolved sources in thelaboratory or on the sky. It is often possible to use geo-metric or physical optics to model the system and prop-agate an image of an unresolved source to the detector.More recently, forward-modeling approaches incorporat-ing measurements of the PSF have been used to modelcomplex electrical effects in detectors (e.g., Donlon et al.2018). In general, astrophysicists have a preference for a r X i v : . [ a s t r o - ph . I M ] F e b Symons et al. using images of unresolved objects to assess the PSFin an image. This method has the advantage of measur-ing the as-built optical system, which not only capturesthe design performance of an optical system, but alsoany effects from scattered light, mechanical structures,or other nonidealities in the telescope or detector.The intrinsic PSF of the optical system and the pix-elization of the displayed image are not immutable quan-tities, but rather independent parameters fixed at eitherthe instrument design, or for the case of time-domainmaps, during the map-making process. This is illus-trated in Figure 1, which shows the effects of changingthe relative sizes of the pixelization and PSF. In prac-tice, the choice of image pixelization depends on thesources under investigation and the desired propertiesof the final image. When possible, the pixel size is cho-sen to be slightly smaller than the intrinsic PSF, but notto “superresolve” it. This is because of both the noisepenalties associated with spreading a fixed flux into anincreased number of pixels, and the diminishing fidelityof superresolved astronomical images. A resolution oftwo pixels per FWHM is often quoted as representativeof Nyquist sampling, with a smaller number of pixels perFWHM indicating undersampling. If the optical PSF issignificantly narrower than the width of a pixel, thenthe PSF is severely undersampled and the precise loca-tion of the source within the pixel may be difficult torecover (Puetter et al. 2005; Robertson 2017). It can beshown that the optimum signal-to-noise ratio (S/N) onunresolved sources is achieved by either concentratingas much of the optical PSF into a single detector ele-ment as possible, or by performing optimal photometrythat uses knowledge of the PSF to weight pixels accord-ing to their flux contributions (e.g., Horne 1986; Naylor1998). The need for accurate photometry of sources isa strong motivation to improve our knowledge of thespatial properties of the PSF, sometimes beyond thatoffered by the native image pixelization.Though the optical PSF may not be resolved by theimage pixelization, information about its shape mustbe retained in any ensemble of images of unresolvedsources that are placed at random locations with re-spect to the pixel grid. In this situation, the images ofthe PSF are sampled by the pixel grid slightly differentlyfor each member of the set, so with a sufficient numberof draws it should be possible to recombine the informa-tion to reconstruct the underlying “unpixelized” opticalPSF. This concept has been successfully used to recon- For our purposes, “unresolved” is the condition θ source (cid:28) θ PSF ,where θ source is the width of the source and θ PSF is the width ofthe PSF. struct superresolved PSFs in the past using a methodcommonly referred to as “stacking” (Dole et al. 2006;B´ethermin et al. 2012; Zemcov et al. 2014). In short,the stacking method is a measurement of the covariancebetween a catalog of sources and their corresponding im-ages in a map. Superresolution methods such as stack-ing can offer a way to recover Nyquist sampling of thePSF (Starck et al. 2002). The stacking method (Cady& Bates 1980) was first used in this context to performstatistical studies on low-resolution, confused images ofthe cosmic infrared background (e.g., Dole et al. 2006;Marsden et al. 2009, and others), but has since beenextended to a wide range of wavelengths and scientifictopics (see, e.g., Viero et al. 2013).In this paper, we study PSF reconstruction and mod-eling using the stacking technique. We demonstratethat, in the limit where the input catalog is composedof point sources, the image resulting from a stackingprocess is that of the intrinsic PSF convolved with thepixelization grid (the effective PSF), and that the pixelgrid can be reliably deconvolved to retrieve the instru-ment’s intrinsic optical PSF. As a case study, we exam-ine this PSF reconstruction algorithm as it may apply tothe Spectro-Photometer for the History of the universe,Epoch of Reionization, and ices Explorer (SPHEREx),an upcoming spectroscopic survey mission designed toimage every 6 . × . pixel over the entire skyfor 0.75 < λ < µ m. Accurate photometry withmethodological uncertainties <
1% is required to helpattain SPHEREx’s ambitious scientific goals, but thePSF is intentionally designed to be underresolved bythe pixel grid with an optical FWHM of 2 . (cid:48)(cid:48)
01 at 0.75 µ m, which is consistent with Korngut et al. (2018) whenmanufacturing tolerance and aberrations are included.This gives a sampling rate of ∼ We define superresolution to mean any resolution higher than animage’s native resolution. econstruction of Severely Undersampled PSFs ∆ φ ( a r c m i n ) A B C ∆ θ (arcmin)3210123 ∆ φ ( a r c m i n ) D ∆ θ (arcmin) E ∆ θ (arcmin) F Figure 1.
Examples of the relationship between pixelization and PSF. The top row shows the effect of changes in the pixelgridding of the input image, shown in panel A, while the bottom row shows the effect of changing the width of the PSF. Inpanel A, the pixelization is matched to the optical PSF, so that the FWHM ∼ θ pix /FWHM ∼ θ pix /FWHM ∼
20 case where the image spatial resolution isheavily gridding-dominated. The method described here takes advantage of the fact that the PSF is sampled in many differentways with respect to the pixel grid to allow reconstruction of the subpixel PSF shape. In the bottom row, we show examples ofFWHM
PSF /θ pix = { , , } in panels D, E, and F, respectively. In these cases, the subpixel PSF can be easily measured frompoint-like sources, and the method described here offers no improvement. of combining multiple exposures to reconstruct the PSFand is less computationally expensive (Seshadri et al.2013). The SPHEREx case illustrates the utility of ourmethod in a situation where the intrinsic PSF is unre-solved in individual images of point sources but can bereconstructed from a large ensemble of measurementscontained within a single exposure. In Section 2, weintroduce the PSF stacking concept in detail, variousmethods of deconvolution necessary for reconstructingthe PSF in complex cases, as well as the method bywhich the PSF reconstruction can be used as a weightfunction in calculating photometry. In Section 3, weapply these methods to example cases including sim-ulated SPHEREx PSFs as well as real data from theLong Range Reconnaissance Imager (LORRI) on NewHorizons (Cheng et al. 2008). In Section 4, we put ourwork in the context of generalized PSF reconstructiontechniques. Table 1 provides a glossary of mathematicalvariables used throughout the text. METHODS2.1.
Stacking
A simple way to understand the relationship betweenthe PSF and the pixel grid function is through the con-volution theorem. The act of signal detection with apixelized system is equivalent to a convolution of thetrue underlying optical PSF, P true , with the pixel re-sponse P grid : P = P true ∗ P grid = F − [ F [ P true ] · F [ P grid ]] , (1)where F represents the Fourier transform and P is thePSF that can be observed in the image itself (Starcket al. 2002; Puetter et al. 2005), sometimes called theeffective PSF (Lauer 1999). Figure 2 illustrates the re-lationship between these quantities.In our formulation, the quantity P true is the combi-nation of any effects that contribute to the overall PSFmeasured in an image. These include the underlyingoptical PSF, how the optical PSF is played around pix-els due to pointing jitter, and any electrical effects like Symons et al.
Table 1.
Important Variables Used in the Text
Variable Description
Image quantities, native pixelization M sky image∆ noise in MP true underlying optical PSF of system P grid pixel grid function P the PSF of sources in M , convolution of P true and P grid Image quantities, super pixelization r pixel grid scaling factor˜ M sky image M scaled up by r ˜ P PSF of sources in M scaled up by r Input source catalog quantities α list of sources in image MF known flux of sources in M ( x, y ) known positions of sources in MN number of sources in M PSF reconstruction quantities T thumbnail cutout of a source centered on its known position P S stacked PSF P RLD deconvolution of P S via RLD P IBP deconvolution of P S via IBP, best estimate for P true I RLD number of RLD iterations I IBP number of IBP iterations
Optimal photometry quantities P P recentered, rescaled version of P IBP used in optimal photometry W P weight function used for optimal photometry F P flux of a source calculated from optimal photometry (cid:104) δF (cid:105) average deviation of flux from its known value σ δF PSF reconstruction figure of merit interpixel capacitance or the “brighter-fatter” effect (Hi-rata & Choi 2020). For a linear detector, the quantity P true incorporates all of these effects as a convolutionover the history of the exposure. It is commonly notedthat the observer only ever has access to P , so it is expe-dient to perform analysis with this quantity. But in thepresence of nonlinear effects that are known to exist inmodern hybridized HgCdTe NIR detectors (Plazas et al.2018; Hirata & Choi 2020), it may in fact be necessaryto deconvolve P grid or other effects from P . There arepotentially many reasons and ways to do this; here, weinvestigate one that uses the fact that the covariancebetween a well-understood catalog of stars and sourcesin the image allows for superresolution reconstruction,which when combined with deconvolution methods canreturn the underlying optical PSF. This is particularlymotivated by the fact that, in the presence of these com- plicating factors, P true is the kernel required to performoptimal photometry.The stacking method takes advantage of the fact thatan “ideal” gridded image of the astronomical sky can berepresented as M j = ∆ j + N j (cid:88) i F i , (2)where M j is the brightness of the image M in pixel j ,∆ j is the noise in that pixel, N j is the number of sourcesfalling into pixel j , and F i is the flux of source i in a list α of N total sources of emission (Marsden et al. 2009;Viero et al. 2013). To make progress with PSF estima-tion, we note that, by the argument above, instrumentswith non-point PSFs contribute flux into more than onepixel. We define P j ( x i , y i ) to be the beam-convolvedand mean-subtracted shape of the PSF centered at some econstruction of Severely Undersampled PSFs Intensity θ θ θ
Power k k k x F F F -1 AD E FB C
Figure 2.
Relationship of the PSF and pixel-gridding function in both real space and Fourier space. Panel A shows ahypothetical PSF, which has a similar size to the pixel-gridding function (dotted line) shown in panel B. The spatial resolutionof the image that results is computed as the convolution of the PSF and a single pixel (solid line in panel B), as shown inpanel C. Flux from a point source is spread across several pixels according to the brightness in the PSF as sampled by the pixelgridding function. As an alternative visualization, panel D shows the power spectrum (absolute square of the Fourier transform)of the PSF, and panel E shows the power spectrum of the pixel-gridding function. By the convolution theorem, the inverseFourier transform of their product, shown as a power spectrum in panel F, is also the final spatial quality of the image of apoint source. These two visualizations are useful for understanding interactions between the PSF and pixel-gridding function. source position ( x i , y i ), and then write Equation 2 as M j = ∆ j + N j (cid:88) i F i P j ( x i , y i ) . (3)This expression accounts for the fact that all sources cancontribute to the intensity of pixel j , as the PSF spreadsflux to neighboring pixels.Our goal is to estimate the shape of the PSF P fromthe measured sky M . Because P has the same amplitudefor each source, we can invert Equation 3 and solve P = N j (cid:88) i M ( x i , y i ) F i , (4)where M ( x i , y i ) is the image centered on the source posi-tion ( x i , y i ), and we assume that the noise obeys (cid:104) ∆ (cid:105) = 0over the sum. Furthermore, as with simple stacking, werequire the source positions to be uncorrelated so thatcontributions from sources in the image but not beingstacked on do not add coherently (Marsden et al. 2009). In this formalism, the superresolution recovery arisesdue to the nature of pixels, which have the property ofaveraging photons. We can always create larger pixelswhich contain more photons, but spread over a largerarea, in such a way that the measured surface bright-ness is conserved. The corollary also holds when goingfrom larger to smaller pixels. The cost of this regriddingoperation is changes in any fixed-amplitude noise: re-gridding larger pixels to smaller ones increases the noiseper pixel, while smaller to larger decreases the noise inthe larger pixel.This averaging property of pixels allows us to writethe relation P j = 1 r r (cid:88) i ˜ P ki , (5)where ˜ P k is the image of the PSF sampled on a finergrid, and the scale factor between the areas of pixels k and j is r . Because the pixel size does not appear in Symons et al.
Equation 4, we can write˜ P = N j (cid:88) i ˜ M ( x i , y i ) F i , (6)where ˜ M is a version of the image on the finer grid andis related to M through M j = (1 /r ) (cid:80) ri ˜ M ki . Any r suchthat the r -times finer image has significantly higher res-olution than the native image resolution will increasethe sampling rate and facilitate the stacking method;increasing r increases the small-scale structure it wouldbe possible to capture in the reconstructed PSF. How-ever, increasing r can also lower the S/N or introducegreater numerical instability due to the larger number ofsubpixels per pixel. The fundamental cost of this super-resolution stacking is that the noise per pixel is increasedover the native pixel resolution. This can easily be ad-dressed by making the list α large to compensate. Asregridding conserves S/N in an area, the total S/N onthe PSF will remain fixed for a fixed noise and numberof sources in the stack.To implement the proposed method in an unbiasedmanner, we require several conditions to be true:1. The source positions are uncorrelated and sourcesare uncrowded, preventing combined or overlap-ping sources from coherently adding in the stack.2. The source list we stack on comprises unresolvedsources so that we are reconstructing an image ofthe PSF.3. The source images do not suffer from detectionartifacts like saturation or nonlinearity.4. (cid:104) ∆ (cid:105) = 0, meaning that stacking over many sourcesaverages down the noise.In practice, we can meet these requirements by makingsuitable choices for the list of sources to stack on andimposing some weak assumptions about the size of thePSF relative to the pixels. Requirements 1 and 2 are metby using a catalog of stars on which to stack. Stars areunresolved by all but very specialized telescopes and (atmiddle and high galactic latitudes where source densityis low) have uncorrelated positions (Zemcov et al. 2014).Requirements 3 and 4 can be met by selecting sourceswith fluxes faint enough to keep in the linear regime ofthe detector and restricting the catalog to a relativelynarrow range of fluxes so the denominator F i in Equa-tion 6 does not strongly overweight noisy sources. How-ever, the noise in the final PSF measurement is a func-tion of the number of sources in the list α , so the trade-off between the flux range and the number of sources depends on the details of both the instrument and thecorresponding survey.In order to approach this problem in a computa-tionally efficient manner, we stack a tractable numberof small “thumbnails” around each source centered on( x i , y i ), which focuses attention on the regions of inter-est and removes the necessity of generating many shiftedversions of M . The size of the thumbnail image is de-termined by including a region large enough that a largefraction of the response from the PSF is included, with-out including unnecessary background noise.With these requirements in mind, the stacking PSFestimation algorithm can be broken into five steps:1. Resample the sky image M into an image ˜ M ona pixel gridding r -times finer. The image will notappear different because a single pixel value from M will fall into multiple pixels in ˜ M , but ˜ M willhave r -times more pixels on a side than M (seeFigure 3).2. For each star in the list α , cut out a thumbnail T centered on the known position of a source ( x, y ).The size of the thumbnail depends on the desiredangular extent of the final PSF estimate. Dueto the nature of resampling, the source’s subpixelcenter will not necessarily be centered in a coarse-grid pixel. The thumbnail is cut such that thesubpixel center will be aligned with the center ofthe stack.3. Suppress the overall constant offset value of thethumbnail. For equally bright sources, subtractingthe mean of the thumbnail image is acceptable,but for sources of varying brightness, an estimateof the sky brightness away from the star is a goodchoice.4. Add the thumbnail centered on the known sourceposition into the stack.5. Repeat for all stars in the catalog. All of the sourcethumbnails are combined by taking the mean ofall corresponding pixels (i.e., the i th pixel in thestack contains the mean of all the i th pixels in eachthumbnail), and this becomes the stacked PSF, P S . This is a superresolution image of the under-lying optical PSF, P true , convolved with the pixelgrid function, P grid . Deconvolution is still neces-sary to return P true (Guillard et al. 2010).Source crowding is a concern when stacking sources.Sources appearing in a thumbnail that have brightnesssimilar to the stacking target will contribute flux tothe stack and broaden the estimate of P . In the limit econstruction of Severely Undersampled PSFs Fine Grid Pixels F i n e G r i d P i x e l s A. High-Resolution Source
Coarse Grid Pixels C o a r s e G r i d P i x e l s B. Downsampled Source
Fine Grid Pixels F i n e G r i d P i x e l s C. Upsampled Source S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] Figure 3.
Regridding of a single source to alter resolution. Panel A shows a single Gaussian point source generated in the finegrid. In panel B, the source resolution is sampled down by a factor of 10 to represent a simulated undersampled image (seeSection 2.4 for details), and in panel C, the resolution is sampled back up by the same factor. Although the source contains thesame number of coarse pixels in panels B and C, the source in panel C contains 100 subpixels in each of these. In panel C, thecenter of the source appears to have shifted from the coarse grid lines due to the subpixel center now being located with moreaccuracy than is possible in the coarse grid. Stacking on each source’s subpixel center makes superresolution PSF reconstructionpossible by allowing the stacked subpixels to generate high-resolution structure.
Table 2.
Source Crowding as a Function of Galactic Latitude
Total Available Sources Isolated Sources Total Sources( (cid:96), b ) with 11 ≤ m AB ≤
15 Sources from Masking in Stack (0 ◦ , 90 ◦ ) 1653 540 808 1348(0 ◦ , 60 ◦ ) 2078 438 1169 1607(0 ◦ , 30 ◦ ) 4710 48 2292 2340(0 ◦ , 15 ◦ ) 6794 5 2185 2190 that source positions are uncorrelated, interloper sourceshave random positions and so act as an extra source ofnoise in P S . However, with a finite number of sourcesthere may be significant sample variance from stack tostack. To mitigate this problem, we mask interlopersources as part of the stacking procedure by measur-ing the distance between the target and any interlopers.Interlopers with center position ≤
11 coarse-grid pixelsfrom the target source’s center are masked by excludingthe thumbnail pixels within a radius of five pixels fromthe interloper position from the sum. If the two sourcesare closer than 8 coarse-grid pixels apart and the tar-get source is not at least an order of magnitude brighterthan the interloper, that source is rejected from the stackentirely. Table 2 demonstrates the effects of our crowd-ing cut on representative 3 . × . (cid:96) , b ) = (0 ◦ , 90 ◦ ) 18% of sources arerejected, while in a crowded field at ( (cid:96) , b ) = (0 ◦ , 15 ◦ )68% of sources are rejected. A benefit of the algorithmdiscussed here is that, even in such crowded fields, weare able to reconstruct a useful estimate for P true (seeSection 3.4). 2.2. Reconstructed PSF Deconvolution
Due to the oversampling inherent to our superresolu-tion stacking method, we require a deconvolution step inorder to remove P grid and return P true from P S (Lauer1999; Starck et al. 2002; Park et al. 2003). This type ofissue has also arisen in Planck data, for which the opticalbeam cannot be recovered without the deconvolution ofthe effect of sampling time (Planck Collaboration et al.2016; Tauber et al. 2019). Guillard et al. (2010) followa similar procedure of superresolution PSF reconstruc-tion followed by deconvolution for the Mid-Infrared In-strument on the James Webb Space Telescope, althoughtheir superresolution method combines multiple imagesfollowed by deconvolution via a maximum a posteriorimethod. We choose the Richardson–Lucy deconvolution(RLD), a common algorithm used in this type of prob-lem (Richardson 1972; Lucy 1974). This algorithm isalso referred to as the expectation-maximization methodand is a form of maximum likelihood estimation (Starcket al. 2002). We have implemented RLD on P S accord-ing to the following prescription. The known blurringfactor P grid (shown in Figure 4) is a matrix where eachpixel contains the fraction of light detected from a point Symons et al. source: P grid ( x, y ) = [ r − ( x − x )][ r − ( y − y )] (7)where r is the upscaling factor, and x and y are thecenter of the stacking area. Starting from P S and P grid ,each iteration of the RLD proceeds as follows: P i +1RLD = P i RLD (cid:18) P S P i RLD ∗ P grid ∗ P refgrid (cid:19) , (8)where the initial P i RLD is simply P S , and i increases onevery iteration for I RLD total iterations. P refgrid is thereflection of P grid across both axes. The final result is P RLD , the deconvolution of P S and a more accurate re-construction of P true . This process is demonstrated forideal Gaussian sources in Figure 4. However, we findthat there is a limit to the resulting improvement in re-construction quality available from increasing the num-ber of RLD iterations due to RLD’s tendency to amplifyartificial structure after too many iterations (Hanischet al. 1997). This is problematic for applications thatrequire fine features of the PSF, for example, optimallyweighted photometry (see Section 2.3).In order to improve the quality of the reconstructedPSF beyond that which can be achieved by RLD, we im-plement a more advanced method of deconvolution withthe iterative back-projection (IBP) algorithm, in whichthe error or difference between simulated and observedlow-resolution images is iteratively reduced (Irani & Pe-leg 1991). This method is also based on maximum like-lihood estimation but was developed for superresolutionimage reconstruction. The approach is similar to back-projection used in tomography and has previously beenused in reconstruction of blurred or degraded images. Itis similar to RLD in that it offers no unique solution andhas the potential to falsely amplify noise after too manyiterations (Park et al. 2003). This method has been ap-plied in remote sensing and planetary science (Bell et al.2006), and other versions of this method have been usedto analyze images of solar flares (Schwartz et al. 2014)and in the testing of superresolution image reconstruc-tion methods on simulated Euclid data (Castellano et al.2015).After P S is obtained, the IBP iterations proceed asfollows:1. An RLD is performed as previously described,with I RLD iterations. This produces P RLD .2. The entire stacking procedure is repeated using P RLD as P true . P RLD is placed into the fine grid(using the same number of sources as were usedpreviously), and the fine grid is downsampled and upsampled by the previously used scale factors.These sources are stacked as before to produce thenew stacked PSF, P NS .3. An error term δP is calculated as δP = P NS − P S . (9)4. δP is used to produce an error-adjusted stack, P AS .On the first iteration, P AS is created via P AS = P S − δP . (10)On every subsequent iteration, P AS is P A( i +1)S = P A( i )S − δP . (11)5. The cycle repeats with the RLD of the newlyformed P A( i +1)S (which is resampled and stackedto form P NS ) until a specified completion criterionis reached. The result of the final iteration is P IBP ,which is the updated and more accurate recon-struction.The overall flow of this algorithm is illustrated in Figure5. Standard RLD fails in the presence of noise for sev-eral reasons, including the amplification of noise overmany iterations and, depending on the choice of imagezero point, the presence of negative values. This arisesbecause the RLD algorithm is based on the maximumlikelihood for Poisson statistics, so its solution requirespositive input data. However, stacking requires a choiceabout the image zero point assumed in each thumbnailthat may produce negative values into P S . Though it ispossible to design algorithms that do not have this fea-ture, the most general case is that we have to accountfor mildly negative values in P S .We handle this problem by implementing a more ad-vanced version of the RLD algorithm (Hanisch et al.1997) that adds an initial estimate of the backgroundand read-noise value to every pixel and suppresses thecontribution from pixels with a value less than thedamping factor in each iteration. This damping pre-vents the amplification of noisy pixels that do not con-tain structure related to the PSF. Because our PSF isconcentrated in the center of the stacked image, we usea damping matrix that performs no damping in the cen-ter but damps heavily beyond a circular radius of 13fine-grid pixels from the center, preventing each iter-ation from drastically changing the values of dampedpixels. Beyond a 13 pixel radius, the PSF is largelynoise dominated and does not contribute much signal tothe optimal photometry. As presented below, this ver-sion of the IBP algorithm successfully deconvolves P grid from P S even in the presence of noise. econstruction of Severely Undersampled PSFs Fine Grid Pixels F i n e G r i d P i x e l s P true Fine Grid Pixels F i n e G r i d P i x e l s P S Fine Grid Pixels F i n e G r i d P i x e l s P grid Fine Grid Pixels F i n e G r i d P i x e l s P RLD S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] Figure 4.
Deconvolution of the pixel grid function. In the top left, P true for a Gaussian source is shown, followed by P S (topright) for a stack of 1000 Gaussian sources at a scaling factor of r = 10, P grid (bottom left) that is computed using Equation 7and deconvolved from P S via RLD, and the resulting P RLD (bottom right). P RLD exactly matches P true , further demonstratedin Figure 7. P SA (i-1) P S P SN P S P RLD P RLD P S RLD ( I RLD
Iterations) thumbnails stacked together P SN P SA δ P First iterationAll other iterationsRepeat processFinal iteration P IBP P RLD P RLD P RLD δ P δ P Figure 5.
Steps of the iterative back-projection algorithm. On the first iteration P S is deconvolved using the regular RLDmethod, while on all subsequent iterations P AS is deconvolved. Deconvolution is followed by a repeat of the stacking procedureusing P RLD as P true , creating P NS . The error is calculated to be the difference between P S and P NS , and that error is subtractedfrom either P S (if first iteration) or the previous P AS , creating the new P AS . The cycle then repeats with deconvolution via RLDagain, resulting in P IBP on the final iteration.
Using a Reconstructed PSF for OptimalPhotometry
Optimal photometry (Naylor 1998) is an example ofan application in which detailed knowledge of the PSF isrequired to reach the maximum possible S/N on pointsource fluxes. The best estimate for P true (here P IBP unless explicitly noted) is shifted via interpolation toaccount for the difference between the source’s known coordinates via catalog reference and those same coordi-nates scaled by the regridding factor, r . It is then down-sampled by r to match the resolution of the source. Thisshifted and resampled reconstructed PSF is normalizedand used to give each pixel the weight of its individualflux contribution. This weight is defined as W P = P P (cid:80) ( P P2 ) , (12)0 Symons et al. where P P is the shifted and resampled reconstructedPSF. Each source’s flux is then computed as F P = (cid:88) ( T · W P ) , (13)where T is a thumbnail cutout of each source, and W P isthe previously calculated weight function. For an imagewith many point sources, we define the average deviationof all source fluxes determined via optimal photometryfrom known (catalog) values to be (cid:104) δF (cid:105) = N (cid:80) i [( F P , i − F i ) /F i ] N , (14)where (cid:104) δF (cid:105) is expressed in percent, N is the number ofsources in the image, F is a source’s known flux, and F P is that computed by optimal photometry.2.4. Image Simulation
In order to test and characterize these algorithms, sim-ulated images of the sky are necessary. We simulate im-ages of undersampled point sources consistent with theSPHEREx expectation (Korngut et al. 2018) accordingto the following prescription. First, we generate a gridthat represents the r -times finer grid or pixelized image.The native detector size for SPHEREx will be 2,048 × . × . field of view (FoV)in each of the six main wavelength bands. We select r to be 10 for Gaussian P true , so that our r -times finerimage contains 20,480 × r to be 20for SPHEREx P true , so that the fine-grid image contains40,960 × ∼ r that we use to simulate the high-resolution image doesnot need to have the same value as that used to increaseresolution during stacking. We choose to set them equalfor simplicity, but for real data, the stacking r can beempirically determined. We populate the generated im-age with point sources using one of two methods:1. One thousand point sources are randomly placedover the FoV, and for simplicity, all sources aregiven uniform flux F ν as a dimensionless quantitysuch that the S/N is chosen to be 20. An S/N of20 corresponds to m AB = 17.87. See the Appendixfor more details on the units we use to define fluxand how we calculate S/N.2. Realistic star fields are generated for any desiredsky position with an all-sky catalog of 332 mil-lion sources, derived by selecting stars in the GaiaDR2 (Gaia Collaboration et al. 2016) catalog with close counterparts (within 1 (cid:48)(cid:48) angular separation)in the AllWISE catalog (Wright et al. 2010). Theuse of the g and r p photometry from Gaia in com-bination with the W1, W2, and W3 photometryfrom WISE gives a rough SED for each source fromwhich realistic fluxes for the SPHEREx bands canbe estimated. For the special case of single-fieldtests presented in Section 3, we use the field cen-tered at the north galactic pole (NGP; ( (cid:96) , b ) =(0 ◦ , 90 ◦ )), where star coordinates are uncorre-lated and fields are the least crowded. This “min-imal” field contains ∼ < m AB < P grid is created, the map is convolved with the P true (ei-ther Gaussian or SPHEREx ) under study. This imageis then sampled down by r = 10 or r = 20 as appropriateto the native image resolution. Figure 3 demonstrateshow the regridding process works for a single source. RESULTS
Table 3.
Summary of Simulations
PSF Type Flux Type RLD/IBP Noise Figures Sections
Gaussian Uniform RLD No 7 3.1Gaussian Catalog RLD No 8 3.2Gaussian Catalog IBP No 9 3.2Gaussian Catalog IBP Yes 10 3.3SPHEREx Catalog RLD No 11,12 3.4SPHEREx Catalog IBP No 11,12,13 3.4SPHEREx Catalog IBP Yes 14,15,16,17,18 3.4
We now have all the pieces required to test the de-scribed method against various cases, including nonidealPSF shapes, noise, and crowding, and to assess its over-all performance. In Sections 3.1 – 3.3, we apply the IBPmethod to Gaussian PSFs in various scenarios to developintuition. In Section 3.4, we introduce the SPHEREx P true and assess the effects of noise, field position, andother quantities of interest. Table 3 gives a summary ofthe various tests and the figures and sections where thecorresponding results can be found. Finally, in Section3.5, we apply the IBP PSF reconstruction to data fromthe LORRI instrument on New Horizons and are able toidentify additional complicating factors present in realdata such as pointing instability. The SPHEREx P true is derived from optical simulations per-formed by L3-SSG as part of the SPHEREx Phase A study usingan end-to-end telescope design. econstruction of Severely Undersampled PSFs Gaussian Point Sources with Uniform Flux
We first explore the fundamental properties of thePSF estimation method using a simple Gaussian modelwithout noise. This allows us to perform an idealizedtest of the various stacking and deconvolution methodswithout complicating factors that may introduce theirown sources of error. We construct a simulated imageas specified in Section 2.4 with Gaussian point sourcesand uniform flux. We begin with the stacking proce-dure described in Section 2.1, which is demonstrated forGaussian point sources in Figure 6. To evaluate thequality of P S as a reconstruction of P true , we comparesections through P S and P true in Figure 7. As expected,we find that P S is significantly wider than P true as P S is P true convolved with P grid .We next evaluate P S ’s performance when used as thekernel for the weighted optimal photometry describedin Section 2.3, and find that the measured fluxes using P S have a wide spread and are >
20% larger than theirexpected values from P true . To improve the accuracyof the photometry, we implement I RLD = 10 iterationsof the RLD algorithm. Figure 7 shows the horizontaland diagonal profiles of P RLD , as well as the differencebetween P RLD and P true , which is negligible comparedto the amplitudes of P S and P true .In order to quantify the accuracy of P RLD as a pixel-weighting kernel for photometry, we simulate 50 realiza-tions of noiseless, constant-flux sources with randomizedsource coordinates. We find (cid:104) δF (cid:105) = 0 . ± . Gaussian Point Sources with Catalog-based Flux
To check the effect of a more realistic distribution ofinput fluxes on the RLD reconstruction, we simulate animage with Gaussian point sources and catalog fluxes fora field at the NGP constructed as described in Section2.4, again without noise in the image. Performing thesame stacking and deconvolution procedure yields theresults shown in Figure 8. All fluxes fall within 1% oftheir known values, but a positive bias at the level of ∼ b > ◦ . The (cid:104) δF (cid:105) from this ensemble is 0.410% ± P RLD still being slightly broader than P true at the level of ∼
1: 10 . We conclude that in order to remove the bias, a more accurate reconstruction of P true is necessary.In order to obtain a more accurate reconstruction, weuse the IBP algorithm outlined in Section 2.2. Similar to I RLD , I IBP is likewise determined by minimizing (cid:104) δF (cid:105) ’stotal spread and bias until no further improvements canbe achieved, which we empirically find converges by I IBP = 10. Results from the more advanced algorithm for the50 independent fields with b > ◦ are shown in Figure9 and show significant improvements in both the meanand variance of δF . This is good evidence that P IBP is amore accurate reconstruction of P true than P RLD . Per-forming the same test on sources with uniform flux alsoresults in the previously seen bias in (cid:104) δF (cid:105) being reducedby an order of magnitude, which further demonstratesthat the bias was due to reconstruction quality and notthe type of source fluxes used.3.3. Noise
We simulate noise from the instrument consistent withSPHEREx Band 1 (centered at 0 . µ m; see Dor´e et al.2018 for full specifications) with pixel RMS of 46 nWm − sr − by adding a white noise component to thenative resolution source image. We also model pho-ton noise from the sources themselves as an additionalsource of noise. At this noise amplitude, a source de-tected at 5 σ has m AB = 19 .
4. We restrict the sourcesincluded in P S to the range 10 ≤ m AB ≤ .
1. Thebright limit is determined by sources that would require >
50% correction due to nonlinearity in the SPHERExdetector, and the lower limit by sources with S/N = 1.The latter choice is not motivated by any known limita-tion of the algorithm, but rather by including only thosesources not dominated by noise. In the NGP field, using P IBP constructed from noisy sources as the kernel for op-timal photometry returns photometric fluxes well withinthe desired 1 % benchmark, as illustrated in Figure 10.Performing a more realistic simulation including noise inboth the photometered source fluxes and the PSF stacksources results in photometry that exceeds the 1 % re-quirement, but that is dominated by the noisy fluxesrather than the intrinsic error from the photometry ker-nel. 3.4.
Spatially Structured PSF
In order to provide a concrete test case and assess ourIBP reconstruction of the SPHEREx PSF, we simulatean image matching SPHEREx’s specifications using thesky catalog for the NGP as described in Section 2.4. Weconvolve the point sources with the SPHEREx P true de-rived from optical simulations of the instrument, shownin Figure 11, and optionally add instrument and pho-ton noise. Next, we perform the stacking procedure,2 Symons et al.
Fine Grid Pixels F i n e G r i d P i x e l s Fine Grid Pixels F i n e G r i d P i x e l s
50 Sources
Fine Grid Pixels F i n e G r i d P i x e l s
500 Sources
Fine Grid Pixels F i n e G r i d P i x e l s S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] S c a l e d F l u x [ d i m e n s i o n l e ss ] Figure 6.
Demonstration of stacking an increasing number of sources. On the top left, five source thumbnails are combinedtogether via their mean. Each new thumbnail is centered on the corresponding source’s known coordinates. The number ofsources increases on the top right to 50 sources, on the bottom left to 500 sources, and on the bottom right to 1000 sources.Although each individual thumbnail bears little resemblance to P true , the more thumbnails that are stacked together, the moreclosely the stack begins to match P true . pixel grid deconvolution, and optimal photometry as de-scribed in Section 2.First, we perform a noiseless test. A comparison of P true to the resulting P S , P RLD , and P IBP calculated inthe absence of noise is displayed in Figure 11. As ex-pected from the Gaussian simulations, P S is not similarto P true , but P RLD with 200 I RLD offers a significant im-provement in reconstruction quality. Using the full IBPalgorithm allows an even more accurate reconstructionof P true . A comparison of the squared Fourier represen-tations of P true and P IBP is shown in Figure 12, alongwith the corresponding differences between P true , P RLD ,and P IBP . Our choice of r dictates the frequency of in-formation we are able to recover in the reconstruction.We find that using r = 10 for the SPHEREx P true doesnot return enough high-frequency information in the re-construction, so we use r = 20 to obtain accurate results.This is an empirical choice based on reconstruction qual-ity as measured by the output photometric accuracy.Applying P IBP as the photometry kernel results in no detectable bias within uncertainties, as demonstrated inFigure 13. With this procedure, all measured fluxes arewithin 1% of their known values, and 50 trials of fieldswith randomized galactic longitude at b > ◦ produce (cid:104) δF (cid:105) consistent with zero.Next, we test the reconstruction of the SPHEREx PSFin the presence of noise in the stacked sources. We findthat setting pixel values beyond an exclusion radius of 31fine-grid pixels to zero after both the initial RLD and theIBP sequence yields the best results in this case. Thisis because noisy pixels with little or no PSF signal canamplify noise artifacts. The damping radius controls theweight given to pixels between the damping and exclu-sion radii during the RLD, while the exclusion radiusprovides a hard cutoff for pixels in the P S image thathave no significant effect on the PSF. Radii for dampingand exclusion are determined through minimization ofthe average flux deviation derived from using P IBP asa kernel for optimally weighted photometry, (cid:104) δF (cid:105) , andits bias from the known input F . Because the average econstruction of Severely Undersampled PSFs Horizontal Cut (Pixels) N o r m a li z e d I n t e n s i t y Diagonal Cut (Pixels) P true P RLD P S P RLD P true Figure 7.
Comparison of P RLD to P S and P true for a stack of Gaussian point sources. The left panel shows a horizontal cutcomparison of the various PSFs and the difference between the input PSF, P true , and the output PSF, P RLD . The horizontalcut (left panel) is a horizontal section with respect to the pixel grid through the center of the PSF, while the diagonal cut (rightpanel) is a section from the top-left to the bottom-right corner. While P S is much broader than P true , the deconvolved P RLD very closely matches the width and shape of the optical PSF.
Figure 8.
The left panel shows the expected vs. measured flux for noiseless simulations of Gaussian PSFs in the field located atthe NGP. Gray lines mark our photometric benchmark within which all fluxes have less than a 1% difference from their knownvalues. Photometry is only performed on sources with m AB >
16 as this marks the (approximate) upper limit of unresolvedgalaxy brightness for SPHEREx. Though the variation is within the requirement, a clear positive bias of the output flux can beseen. The right panel shows the cumulative distribution function (CDF) of (cid:104) δF (cid:105) for 50 trials of fields with b > ◦ . The meandeviation, 0 . ± . P true . deviation from the expected flux is desired to be zero,any significant trend from zero indicates a bias beingintroduced during reconstruction or photometry. Thevalue and bias of (cid:104) δF (cid:105) are minimized until the numberof iterations, I IBP , and the exclusion and damping radiiboth converge. The radii and number of iterations haveconverged when no further improvements occur. Fig- ure 14 demonstrates the deconvolution process for theSPHEREx PSF in the presence of noise.3.4.1.
Evaluating Reconstruction Quality
In order to evaluate the quality of the reconstruction,as well as its effectiveness when used as a weight kernelfor optimal photometry, we develop a figure of merit(FOM) based on (cid:104) δF (cid:105) . Our FOM is the RMS of δF :4 Symons et al.
Figure 9.
The same test as shown in Figure 8 is repeated with the addition of 10 I IBP . On the left, the comparison of thedifference between measured and expected flux to expected flux for a single field at the NGP reveals a much tighter correlationwith bias at the 10 − level. On the right, the CDF of (cid:104) δF (cid:105) from 50 trials with varying galactic longitude at b > ◦ has a meanof 0 . ± . Figure 10.
For a single trial using a Gaussian P true and catalog-based sources, varying photometric accuracy results dependingon whether we add noise to the sources, the kernel, or both. Panel A shows the photometry using a noisy P IBP as the kernelfor noisy sources. Results are similar in panel B, which shows flux values obtained using the noiseless P true applied to noisysources. For both panels A and B, instrumental and photon noise dominate the photometric results. Using the noisy P IBP asthe kernel for photometry of noiseless sources (panel C) gives all flux values within the desired 1% requirement, demonstratingthat our reconstruction methods are not introducing photometric error beyond the SPHEREx requirement. σ δF = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) N (cid:80) i [( δF i ) ] N (15)where δF i is the difference between measured and ex-pected flux for any individual source, and N is the totalnumber of sources. This gives a measure of the total dis-persion in δF such that σ δF < P S . Bright sources allow mea-surement of P S in the coarse gridding with high fidelity,but as there are relatively few sources, this comes at the cost of the spatial resolution of the reconstruction.Faint sources are numerous, but the effect of noise islarger in the stack, and this has a cascading effect onthe fidelity of P IBP . In order to optimize the range offluxes to use, we calculate P IBP over narrow magnituderanges and from these calculate σ δF to determine whereit is optimal. Figure 15 demonstrates how the average σ δF of 50 trials with randomized galactic longitude for b > ◦ changes when selecting only sources within asingle AB magnitude bin for use in P S . We find that σ δF is minimized for 11 ≤ m AB ≤
15, which are thestatistically optimal source magnitudes for optimal PSFreconstruction at the SPHEREx noise amplitude. econstruction of Severely Undersampled PSFs Fine Grid Pixels F i n e G r i d P i x e l s P true Fine Grid Pixels F i n e G r i d P i x e l s P S Fine Grid Pixels F i n e G r i d P i x e l s P RLD
Fine Grid Pixels F i n e G r i d P i x e l s P IBP I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] Figure 11.
Images of the SPHEREx P true (top left), P S (top right), P RLD after 200 RLD iterations (bottom left), and P IBP after 300 RLD iterations combined with the adaptive IBP iterations (bottom right) on 0 . (cid:48)(cid:48)
31 pixel − gridding. Because the IBPoffers a more accurate reconstruction of P true , it also provides more accurate photometry when used as pixel weighting. Next, we test how using a much smaller subset of thetotal available sources affects the reconstruction qual-ity. We restrict the number of sources used in P S byan increasing percentage via a representative sample forthe optimal range of 11 ≤ m AB ≤
15. Figure 16 showshow the average σ δF varies with this restriction for aseries of 50 trials of fields with varying galactic longi-tude and b > ◦ at each fraction of the total num-ber of sources. As the fraction of sources being usedincreases, σ δF continues to decrease, indicating betterreconstruction quality. Even restricting to only 6% ofthe available sources yields σ δF < P S appropriate to a region wherethe gradient of the PSF properties will be smaller.3.4.2. Overall Photometric Accuracy
To assess the overall reliability of the IBP algorithm,in Figure 17 we compare all six SPHEREx wavelengthbands (centered at 0.930, 1.375, 2.030, 3.050, 4.065, and4.730 µ m) using the 11 ≤ m AB ≤
15 range for P S . P true and the noisy P IBP are shown for each band, along with a cumulative distribution function (CDF) of theresulting (cid:104) δF (cid:105) for the noisy P IBP applied to noiselesssources for 50 trials with randomized galactic longitudeat b > ◦ . We use catalog-based fluxes from the Gaiacatalog for all bands, resulting in a high level of accu-racy with mean percent error consistent with zero foreach band. This demonstrates the algorithm’s robust-ness in successful PSF reconstruction of complex PSFswith different types of structure from images with vary-ing levels of noise.In previous tests, we have performed optimal photom-etry on intrinsically noiseless sources to assess the ac-curacy of the IBP reconstruction algorithm. It is alsoinformative to investigate at which source fluxes the in-strumental and shot noise present in real sources woulddominate over the error arising from misestimation of P true . In the top row of Figure 18, we show tests usingcombinations of noise in the stack and noise in the pho-tometered sources (or lack thereof) for the NGP, com-pared to using a noiseless P true as the photometry ker-nel. We continue to restrict P S to 11 ≤ m AB ≤
15. Asexpected, only when the noisy P IBP is used for photome-try on sources with perfectly known flux is the resultingphotometry within 1% variation at all F . Photometryon noisy sources is dominated by the noise in the sourceflux measurement. Noise in the PSF reconstruction it-6 Symons et al.
20 10 0 10 20 k x [1/arcsec] k y [ / a r c s e c ] F [ P true ]
20 10 0 10 20 k x [1/arcsec] k y [ / a r c s e c ] F [ P IBP ]
20 10 0 10 20 k x [1/arcsec] k y [ / a r c s e c ] F [ P true P RLD ]
20 10 0 10 20 k x [1/arcsec] k y [ / a r c s e c ] F [ P true P IBP ] [ d B ] [ d B ] [ d B ] [ d B ] Figure 12.
A comparison between the squared Fourier transform of the SPHEREx P true (top left), P IBP (top right), differencebetween P true and P RLD (bottom left), and difference between P true and P IBP (bottom right). We would expect a diffraction-limited system to go to zero response near 4 πD/λ ∼
14 arcsec − ; the measured value of 11 . − is consistent with thebeam being slightly larger than diffraction limited. The IBP method exhibits residuals at the 10 − level, which is significantlyless than the traditional RLD method that shows residuals at the > − level. Figure 13. (Left) The difference between measured and expected flux as a function of expected value, with gray lines showingthe 1% requirement for a single trial using the noiseless SPHEREx P true and IBP algorithm. All fluxes lie within the 1%boundaries, and we detect no overall flux bias. (Right) (cid:104) δF (cid:105) for 50 trials of fields at b > ◦ and randomized galactic longitude,which has a mean of − . ± . self appears to have a small effect compared with theintrinsic scatter of F from random noise. Additionally,the PSF reconstruction is not correlated with details ofthe noise realization; the P IBP derived from one noise realization, when used for photometry on sources fromdifferent noise realizations, returns an unbiased (cid:104) δF (cid:105) .To understand the relative effect of noise in the stackversus numerical inaccuracies in the algorithm, we per- econstruction of Severely Undersampled PSFs Fine Grid Pixels F i n e G r i d P i x e l s P true Fine Grid Pixels F i n e G r i d P i x e l s P S Fine Grid Pixels F i n e G r i d P i x e l s P RLD
Fine Grid Pixels F i n e G r i d P i x e l s P IBP I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] Figure 14.
The deconvolution process in the presence of noise for the 0.93 µ m SPHEREx wavelength band. In the upper-leftpanel we show P true for comparison with the noisy P S on the top right, the initial P RLD with the use of damping and backgroundnoise constant on the bottom left, and the final P IBP after the IBP process on the bottom right. Due to the large number ofhigh-S/N sources in the stack, instrument noise has minimal effect on reconstruction quality. form 50 trials with fields of varying galactic longitudeat b > ◦ for the same set of noise cases. We againfind that the variation in the output flux is dominatedby the intrinsic noise rather than the error in the IBP,demonstrated in the bottom row of Figure 18. Resultsare very similar for noisy sources with a noiseless P true and noiseless sources with a noisy P IBP , as both havemean (cid:104) δF (cid:105) consistent with zero and similar distributionwidths. Again, the PSF kernel reconstruction is not in-troducing photometric error beyond that determined bythe noise on the sources.3.4.3. Source Crowding
To understand the effects of source crowding (dis-cussed in Section 2.1), we vary galactic coordinates ina single trial for each of b = 15 ◦ , 30 ◦ , 60 ◦ , and (cid:96) = 0 ◦ ,90 ◦ , and 180 ◦ with P S restricted to 11 ≤ m AB ≤ P true for the 0.93 µ m wave-length band with catalog-based fluxes and noisy P IBP with noiseless sources during photometry. Table 4 givesa comparison of (cid:104) δF (cid:105) and σ δF for each set of coordinates.All fields have σ δF < (cid:104) δF (cid:105) and no obvi-ous increase in σ δF . This indicates that accurate PSF re-construction and photometry can be achieved with thismethod, even for highly crowded fields. Table 4.
Photometric Results as a Function of Sky Position (cid:96) = 180 ◦ (cid:96) = 90 ◦ (cid:96) = 0 ◦ (cid:104) δF (cid:105) σ δF (cid:104) δF (cid:105) σ δF (cid:104) δF (cid:105) σ δF b = 60 ◦ b = 30 ◦ –0.046% 0.201% 0.004% 0.168% –0.026% 0.207% b = 15 ◦ –0.078% 0.205% –0.066% 0.172% –0.090% 0.210% LORRI PSF Estimation
To this point, we have been working solely with sim-ulated PSFs where we know the input and the noise isidealized. To test this algorithm against real data, weoperate on an image taken by the LORRI instrument onNew Horizons (Cheng et al. 2008; Zemcov et al. 2017).This image was acquired by the CCD chip on the LORRIinstrument and, in LORRI’s 4 × ×
256 pixels at 4 . (cid:48)(cid:48) − , as shown in Fig-ure 19. In the PSF reconstruction, we use r = 20 and I IBP = 10. The reconstructed PSF is shown in Figure19. We find that the FWHM of the minor axis of P IBP is 3 . (cid:48)(cid:48)
73, while the FWHM of the major axis is 5 . (cid:48)(cid:48) . (cid:48)(cid:48) Symons et al. m AB F [ % ] S/N N u m b e r o f S o u r c e s Figure 15.
Comparison of σ δF when restricting the sourcesused in P S to 1 mag wide bins (lower axis) centered betweeninteger m AB values. The blue points and line indicate theaverage σ δF (left axis) from 50 trials of each bin at varyinggalactic longitude and b > ◦ . The black dashed line showsthe desired limit of σ δF = 1, and the green dashed line in-dicates the average σ δF from 50 trials of a 11 ≤ m AB ≤ σ error. This regionprovides the lowest σ δF compared to individual magnitudebins. The gray line shows the number of sources used in P S in any given bin (right axis), and the upper axis shows thecorresponding S/N for each magnitude bin. σ δF always de-creases with an increasing number of sources or higher S/Nand is minimized in the zone that balances both. Fraction of Sources [%] F [ % ] Figure 16.
Comparison of σ δF correlating to the fractionof the total number of available sources being used in P S .The blue region corresponds to the 1 σ error from 50 trialsof randomized galactic longitude at b > ◦ . Increasing thenumber of sources used in P S results in lower σ δF and thusbetter reconstruction quality, while using as few as 6% ofsources still allows for σ δF <
1, marked by the black dashedline.
What can account for this difference? Noble et al.(2009) demonstrate that the PSF width is a strong func-tion of the integration time of the instrument, and theNew Horizons spacecraft is known to exhibit pointingdrift at the arcsec s − level. Performing the deconvo-lution on a set of several hundred 10 s LORRI expo-sures, we find the PSF is often extended with eccentric-ity (cid:15) > . (cid:48)(cid:48) . Thisis consistent with the expected drift in the spacecraft’spointing in relative control mode of ± (cid:48)(cid:48) per exposure(Conard et al. 2017). In principle, if we understood thepointing history of the spacecraft, we could deconvolvethis component out to isolate the underlying optical PSFusing the same methods described above; this will be leftto future work. DISCUSSION4.1.
Comparison to Similar Methods
The method of PSF reconstruction described here isrobust to complicating factors such as severe undersam-pling, complex PSFs, noise, crowded fields, or a lim-ited number of available point sources. Starck et al.2002 and Puetter et al. 2005 provide thorough reviewsof deconvolution methods commonly used in astronomy,which are numerous. New methods are still being ac-tively developed for applications in astronomy (see, e.g.,Sureau et al. 2020). In comparison to some algorithms,a major advantage of this method is that it works on aper-exposure basis. Many other superresolution meth-ods that have been used in astronomy (e.g., Orieux et al.2012; Li et al. 2018; Guo et al. 2019) rely on multiple ex-posures of the same field for information reconstruction.As an example, data from the Hubble Space Telescope(HST) prompted the development and application of anumber of techniques driven by undersampling in theWide Field and Planetary Camera 2 (WFPC2; Ander-son & King 2000) and Wide Field Camera 3 (WFC3;Anderson 2016). These methods rely on dither-basedsolutions in which multiple images of the same field withsome small shift in the position of the detector are re-combined to generate a higher-resolution image of thesame field (Lauer 1999; Fruchter & Hook 2002). Ourmethod does not rely on similar constraints.The SPHEREx application investigated here requiressuperresolution knowledge of the PSF to optimallyweight pixels for photometry and, as a result, requiresthe effect of the pixelization to be modeled on a per-source basis. Some reconstruction methods return thePSF convolved with the pixel-gridding function (the ef-fective PSF or the PSF observed in the image) insteadof the underlying optical PSF, such as the method in-troduced by Anderson & King (2000). While the effec- econstruction of Severely Undersampled PSFs Fine Grid Pixels F i n e G r i d P i x e l s B a n d True PSF
Fine Grid Pixels F i n e G r i d P i x e l s Reconstructed PSF -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) PSF-Weighted Photometry
Fine Grid Pixels F i n e G r i d P i x e l s B a n d Fine Grid Pixels F i n e G r i d P i x e l s -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] Fine Grid Pixels F i n e G r i d P i x e l s B a n d Fine Grid Pixels F i n e G r i d P i x e l s -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) Fine Grid Pixels F i n e G r i d P i x e l s B a n d Fine Grid Pixels F i n e G r i d P i x e l s -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] Fine Grid Pixels F i n e G r i d P i x e l s B a n d Fine Grid Pixels F i n e G r i d P i x e l s -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) Fine Grid Pixels F i n e G r i d P i x e l s B a n d Fine Grid Pixels F i n e G r i d P i x e l s -0.4 -0.2 0.0 0.2 0.4 F from 50 Trials p ( F X ) I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] I n t e n s i t y I [ n W m s r ] Figure 17.
Comparison of the six SPHEREx wavelength bands, with each band’s P true in the first column, final reconstructed P IBP in the second column, and CDF of (cid:104) δF (cid:105) from 50 trials of randomized galactic longitude at b > ◦ in the third column.Mean (cid:104) δF (cid:105) for all trials is marked with a dotted line, with a gray region showing the standard error of the mean. The mean (cid:104) δF (cid:105) is consistent with zero for each band. Symons et al.
Figure 18.
The top row shows single trials at the NGP using the SPHEREx 0.93 µ m P true with catalog-based flux values andinstrument noise. Panels A and B show photometry resulting from noisy sources with a noisy P IBP and noisy sources with anoiseless P true , respectively. In both cases, flux values display similar scatter, indicating that reconstruction quality does notnegatively affect photometric results. Panel C shows the results from noiseless sources and a noisy P IBP , all of which are within1% boundaries. To quantify the bias in these cases, the bottom row shows CDFs of (cid:104) δF (cid:105) corresponding to panels A – C for 50trials with randomized field galactic longitude and b > ◦ , where the dotted line indicates the mean (cid:104) δF (cid:105) of all trials, and thegray region shows the standard error of the mean. Panel D shows a mean (cid:104) δF (cid:105) of 0 . ± . (cid:104) δF (cid:105) of − . ± . (cid:104) δF (cid:105) of 0 . ± . Figure 19.
On the left is an example LORRI image with 4. (cid:48)(cid:48) − . On the right is a comparison of the P S and P IBP measured from this image, shown in linear scale in the top row and logarithmic scale in the bottom row. The PSF measuredin this image is significantly larger than the laboratory-determined PSF (Cheng et al. 2008), which is likely due to the pointingstability of the spacecraft (Noble et al. 2009). econstruction of Severely Undersampled PSFs ∼ Future Improvements
To summarize the primary result of this work, throughsimulations of the stacking method in all six SPHERExwavelength bands including realistic noise, source cat-alogs from Gaia+AllWISE star catalogs, and realisticbeam shapes, we find IBP-derived kernels allow pho-tometry with accuracy to ∼ ≤ m AB ≤
15 generate the best estimate of the un-derlying optical PSF across all six SPHEREx bands. AtSPHEREx’s noise level, the population of sources in thisrange balance the need for large S/N in the stack againstthe need for a large number of sources to maintain spa-tial fidelity. We find that stacking on subimages of thefull SPHEREx FoV does not significantly degrade theperformance of the algorithm up to fields with (cid:46)
Software:
Astropy (Astropy Collaboration et al.2013, 2018), Matplotlib (Hunter 2007), NumPy (VanDer Walt et al. 2011), SciPy (Virtanen et al. 2020).2
Symons et al.
REFERENCES
Anderson, J. 2016, Empirical Models for the WFC3/IRPSF, Space Telescope WFC Instrument Science ReportWFC3 2016-12Anderson, J., & King, I. R. 2000, PASP, 112, 1360,doi: 10.1086/316632Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz, B. M.,et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4fAujol, J.-F., Gilboa, G., Chan, T., & Osher, S. 2006,International Journal of Computer Vision, 67, 111,doi: 10.1007/s11263-006-4331-zBell, J. F., Joseph, J., Sohl-Dickstein, J. N., et al. 2006,JGRE, 111, E02S03, doi: 10.1029/2005JE002444Bertin, E. 2011, in Astronomical Society of the PacificConference Series, Vol. 442, Astronomical Data AnalysisSoftware and Systems XX, ed. I. N. Evans,A. Accomazzi, D. J. Mink, & A. H. Rots, 435B´ethermin, M., Le Floc’h, E., Ilbert, O., et al. 2012, A&A,542, A58, doi: 10.1051/0004-6361/201118698Cady, F. M., & Bates, R. H. T. 1980, OptL, 5, 438,doi: 10.1364/OL.5.000438Castellano, M., Ottaviani, D., Fontana, A., et al. 2015, inASPC, Vol. 495, Astronomical Data Analysis Softwareand Systems XXIV (ADASS XXIV), ed. A. R. Taylor &E. Rosolowsky, 257. https://arxiv.org/abs/1501.03999Cheng, A. F., Weaver, H. A., Conard, S. J., et al. 2008,SSRv, 140, 189, doi: 10.1007/s11214-007-9271-6Conard, S. J., Weaver, H. A., N´u˜nez, J. I., et al. 2017, inSociety of Photo-Optical Instrumentation Engineers(SPIE) Conference Series, Vol. 10401, Proc. SPIE,104010W, doi: 10.1117/12.2274351Dole, H., Lagache, G., Puget, J.-L., et al. 2006, A&A, 451,417, doi: 10.1051/0004-6361:20054446Donlon, K., Ninkov, Z., & Baum, S. 2018, PASP, 130,074503, doi: 10.1088/1538-3873/aac261Dor´e, O., Werner, M. W., Ashby, M. L. N., et al. 2018,arXiv e-prints, arXiv:1805.05489.https://arxiv.org/abs/1805.05489Fish, D. A., Brinicombe, A. M., Pike, E. R., & Walker,J. G. 1995, JOSAA, 12, 58,doi: 10.1364/JOSAA.12.000058Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144,doi: 10.1086/338393Gai, M., & Cancelliere, R. 2007, MNRAS, 377, 1337,doi: 10.1111/j.1365-2966.2007.11693.xGaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al.2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272 Guillard, P., Rodet, T., Ronayette, S., et al. 2010, in SpaceTelescopes and Instrumentation 2010: Optical, Infrared,and Millimeter Wave, ed. J. M. O. Jr., M. C. Clampin, &H. A. MacEwen, Vol. 7731, International Society forOptics and Photonics (SPIE), 166 – 178,doi: 10.1117/12.853591Guo, R., Shi, X., & Wang, Z. 2019, JEI, 28, 023032,doi: 10.1117/1.JEI.28.2.023032Hanisch, R. J., White, R. L., & Gilliland, R. L. 1997, inDeconvolution of Images and Spectra: Second Edition,ed. P. A. Jansson (San Diego, CA: Academic Press, Inc.),310–360Hirata, C. M., & Choi, A. 2020, PASP, 132, 014501,doi: 10.1088/1538-3873/ab44f7Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801Hunter, J. D. 2007, CSE, 9, 90, doi: 10.1109/MCSE.2007.55Irani, M., & Peleg, S. 1991, CVGIP: Graphical Models andImage Processing, 53, 231,doi: 10.1016/1049-9652(91)90045-LKorngut, P. M., Bock, J. J., Akeson, R., et al. 2018, inSociety of Photo-Optical Instrumentation Engineers(SPIE) Conference Series, Vol. 10698, Proc. SPIE,106981U, doi: 10.1117/12.2312860Lauer, T. R. 1999, PASP, 111, 1434, doi: 10.1086/316460Li, Z., Peng, Q., Bhanu, B., Zhang, Q., & He, H. 2018,Ap&SS, 363, 92, doi: 10.1007/s10509-018-3315-0Lucy, L. B. 1974, AJ, 79, 745, doi: 10.1086/111605Marsden, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ,707, 1729, doi: 10.1088/0004-637X/707/2/1729Naylor, T. 1998, MNRAS, 296, 339,doi: 10.1046/j.1365-8711.1998.01314.xNoble, M. W., Conard, S. J., Weaver, H. A., Hayes, J. R.,& Cheng, A. F. 2009, in Society of Photo-OpticalInstrumentation Engineers (SPIE) Conference Series,Vol. 7441, Proc. SPIE, 74410Y, doi: 10.1117/12.826484Orieux, F., Giovannelli, J. F., Rodet, T., et al. 2012, A&A,539, A38, doi: 10.1051/0004-6361/201116817Park, S., Park, M., & Kang, M. 2003, ISPM, 20, 21 ,doi: 10.1109/MSP.2003.1203207Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016,A&A, 594, A7, doi: 10.1051/0004-6361/201525844Plazas, A. A., Shapiro, C., Smith, R., Huff, E., & Rhodes, J.2018, PASP, 130, 065004, doi: 10.1088/1538-3873/aab820Puetter, R., Gosnell, T., & Yahil, A. 2005, ARA&A, 43,139, doi: 10.1146/annurev.astro.43.112904.104850Richardson, W. H. 1972, JOSA, 62, 55,doi: 10.1364/JOSA.62.000055Robertson, J. G. 2017, PASA, 34, e035,doi: 10.1017/pasa.2017.29 econstruction of Severely Undersampled PSFs Rowe, B., Hirata, C., & Rhodes, J. 2011, ApJ, 741, 46,doi: 10.1088/0004-637X/741/1/46Schmitz, M. A., Starck, J. L., Ngole Mboula, F., et al.2020, A&A, 636, A78, doi: 10.1051/0004-6361/201936094Schwartz, R. A., Torre, G., & Piana, M. 2014, arXive-prints, arXiv:1407.7343.https://arxiv.org/abs/1407.7343Seshadri, S., Shapiro, C., Goodsall, T., et al. 2013, PASP,125, 1065, doi: 10.1086/673318Shi, X., Guo, R., Zhu, Y., & Wang, Z. 2017, Journal ofSystems Engineering and Electronics, 28, 1236,doi: 10.21629/JSEE.2017.06.21Starck, J. L., Pantin, E., & Murtagh, F. 2002, PASP, 114,1051, doi: 10.1086/342606 Sureau, F., Lechat, A., & Starck, J.-L. 2020, A&A, 641,A67, doi: 10.1051/0004-6361/201937039Tauber, J. A., Nielsen, P. H., Mart´ın-Polegre, A., et al.2019, A&A, 622, A55, doi: 10.1051/0004-6361/201833150Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011,CSE, 13, 22, doi: 10.1109/MCSE.2011.37Viero, M. P., Moncelsi, L., Quadri, R. F., et al. 2013, ApJ,779, 32, doi: 10.1088/0004-637X/779/1/32Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020,Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al.2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868Zemcov, M., Immel, P., Nguyen, C., et al. 2017, NatureCommunications, 8, 15003, doi: 10.1038/ncomms15003Zemcov, M., Smidt, J., Arai, T., et al. 2014, Science, 346,732, doi: 10.1126/science.1258168
APPENDIX A. POINT-SOURCE FLUX AND UNITSFor point sources that are assigned some magnitude and corresponding specific flux F λ (measured in units of powerper unit area per unit wavelength such as nW m -2 µ m -1 ), we calculate specific intensity (generally measured in unitsof power per unit area per unit solid angle per unit wavelength) I λ as I λ = F λ Ω beam [nW m − sr − µ m − ] (A1)where Ω beam is the beam area on the sky for SPHEREx. We then multiply I λ by a specific wavelength to get theintensity or diffuse surface brightness λI λ , where λI λ = λF λ Ω beam [nW m − sr − ] . (A2)The quantity λI λ is equivalent to the quantity νI ν , where I ν is measured in units of power per unit area per unit solidangle per unit frequency such as nW m -2 sr -1 Hz -1 . With the proper conversion, I λ and I ν can be used interchangeably.When we perform photometry of a point source, we are taking the integral (cid:90) beam I ν d Ω = F ν [Jy] (A3)to get some flux F ν . For our point sources with catalog-based flux values, we choose to convert the photometric fluxinto intensity in units natural to SPHEREx, nW m -2 sr -1 , using equation A2. Given some source with intensity λI λ and a pixel RMS σ , the S/N can be defined as S / N = λI λ σ (A4)where λI λ and σ have been integrated over a beam with area on the sky Ω beam . For our simulated point sources withuniform flux, we assign λI λλ