Astronomical Image Processing with Array Detectors
ABSTRACT
We address the question of astronomical image processing from data obtained with array detectors. We define and analyze the cases of evenly, regularly, and irregularly sampled maps for idealized (i.e., infinite) and realistic (i.e., finite) detectors. We concentrate on the effect of interpolation on the maps and the choice of the kernel used to accomplish this task. We show how the normalization intrinsic to the interpolation process must be carefully accounted for when dealing with irregularly sampled grids. We also analyze the effect of missing or dead pixels in the array and their consequences for the Nyquist sampling criterion.
Received 2007 April 23; accepted 2007 June 25; published 2007 August 20
1. INTRODUCTION
The creation of smooth two‐dimensional maps from a series of samples measured at discrete points is a common problem in astronomical image processing. The goal is to create a smooth map with the best possible spatial resolution, given a set of data sampled in two dimensions. The solution is complicated by the fact that the data are often not sampled in a regular way, even if the detector layout is regular. For example, telescopes may be scanned or dithered to map areas larger than the array, some instruments are unable to follow objects as they rotate on the sky, or an array itself may contain flaws (i.e., missing or dead pixels). The resulting two‐dimensional sample pattern can often appear quite irregular.
Although the layout of most modern detector arrays (e.g., CCDs) can reasonably be approximated as generating evenly sampled grids (ESGs) extending to infinity in all directions, observations with these arrays will typically include a series of array translations and rotations. In addition, one may want to combine multiple images of the same piece of sky in order to increase the signal‐to‐noise ratio (S/N). This requires that the images be registered so that the same area of the sky is being observed in each image. That is, any relative translation and rotation of the array positions with respect to the sky must be taken into account when combining the images. Unless the translation and rotation operations are such that every pixel lies in a location previously occupied by another pixel, the resulting sample pattern is no longer an ESG.
The processing of data from any case other than an ESG requires performing an interpolation. Some of these cases have been discussed by other authors (e.g., Granrath & Lersch 1998). The interpolation (or smoothing) necessarily has an effect on the spatial resolution of the resulting map. For the case of any sampling pattern (ESG or otherwise), we wish to address the following two questions: (1) how does one choose an optimal kernel shape and size for the interpolation function, and (2) how does this kernel choice affect the spatial resolution of the resulting map? We concentrate on the case of a Gaussian kernel.
As mentioned above, the construction of maps from non‐ESG sampled data is generally done through an interpolation of the data using a smoothing kernel (see § 4 and Lombardi & Schneider 2001). In this paper, we begin by reviewing the solution for ESGs using a technique based on Fourier transforms (§ 3) and extend this technique in § 4 to regularly spaced grids (RSGs), which are composed of relatively translated ESGs. In § 5, we use the tools developed for studying RSGs and ESGs to analyze irregularly sampled grids (ISGs) and explore the effects of missing samples in a map.
Throughout this paper, we present examples that reference the bolometer array used with SHARP, the Submillimeter High Angular Resolution Polarimeter. SHARP is used in conjunction with the Submillimeter High Angular Resolution Camera II (SHARC‐II), which is deployed at the Caltech Submillimeter Observatory (Dowell et al. 2003). For this, the 
SHARC‐II detector array is optically split into two 
subarrays (and a section of 
unused pixels), which image two orthogonal linear polarization components of radiation (Novak et al. 2004). Although we concentrate on SHARP maps, the results are applicable to any detector array or sampling pattern.
2. MATHEMATICAL DEFINITIONS
Before embarking on the analysis of the evenly sampled grid (ESG), we first introduce a set of definitions and functions that are central to the development of the subsequent sections. Given a function 
, which is dependent on position 
, (
and 
are the usual Cartesian unit basis vectors), we define the Fourier transform pair1 
where 
is the spatial frequency vector. The digitization of a signal will invariably introduce trains of Dirac distributions. For example, a two‐dimensional Dirac train of periods 
and 
along 
and 
, respectively, is defined such that 
with the Fourier transform relation (see the Appendix) 
where 
Another useful distribution is the flat‐top window of length 
(in the 
direction in this case), which we denote by 
and the corresponding Fourier transform pair 
3. THE EVENLY SAMPLED GRID
The detection of a signal 
from an astronomical source is inevitably achieved through a series of transformations. Mathematically speaking, the signal is first convolved with the telescope transfer function 
such that 
where the symbol “
” stands for a convolution, while the measured signal 
is a sampled, pixel‐integrated version of 
. For an ESG, the sampling is done in an even manner with a Dirac train as defined in equations (3) and (4). More precisely, for rectangular pixels of widths 
and 
, we write 
with 
The convolution 
stands for what is commonly described as the point‐spread function (PSF). Using equations (2) and (4) and the properties of the Fourier transform for products and convolutions of functions, we find that 
with 
and 
Equation (10) is only valid for the idealized case of an infinite array. In reality, this relation should be multiplied by an aperture function of appropriate width and shape. Although we take this restriction into account when analyzing the effect of missing pixels in § 5.2, we will for the moment simplify our analysis by assuming that the array is sufficiently large so that equations (10) and (14) are suitable approximations.
3.1. Interpolation
The ESG studied in the previous section is the simplest representation that can be given for a sampled set of data. As we see in later sections, we will always seek to transform more complicated forms of data grids (i.e., not evenly sampled ones) into ESGs to facilitate analysis; this will invariably require the interpolation of sampled quantities from different locations. In addition, one might inquire about quantities at positions where there are no samples. For example, questions such as “What is the intensity at position A, where there is no sample, and how does it compare to the flux at positions B, C, and D on this map?” are common when analyzing astronomical images. It is therefore often necessary to generate a new interpolated map from the data set expressed through equation (10).
Given a weighting function 
, any value 
to be assigned to an interpolated point can be expressed as 
where 
is the value associated with the ith of the n data points used for the interpolation. The quantity 
is the normalization factor, which is a function of the position of interpolation. The generation of an interpolated map is equivalent to the convolution of the initial data set with the weighting function, followed by the normalization and resampling of the data. This can be ascertained through a comparison of equation (18) with 
where 
is defined as in equation (5) and 
is the displacement vector specifying the position of the interpolated grid 
in the relation to the initial grid. It is important to realize that because of the evenness in the sampling distribution of the original map 
, the normalization factor 
will be periodic in character, with the same periods (i.e., 
and 
) as the original sampling Dirac train.2 As a consequence, it will take a common value for all interpolated points similarly located within a one‐period segment anywhere on the grid (see the Appendix). More precisely, data resulting from interpolations at points 
and 
, for any integer i and k when 
is defined as in equation (5), will have the same normalization factor. Therefore, when the resampling is done using the same spatial sampling rate as for the original grid (as is the case in eq. [20]), we can write Fourier transform of 
as 
where c is the constant value associated with 
for this particular resampling process. Equation (22) contains multiple copies of 
, one for each pair of m and n. If the Nyquist sampling criterion is satisfied (see the Appendix), then the high‐frequency copies may be removed with negligible aliasing by choosing the weighting function such that 
when 
(when 
or 
). Equation (22) then simplifies to 
and 
The only difference between 
and 
(see eq. [10]), besides the overall scaling factor and translation, is the presence of the convolution by 
for the former. It is therefore apparent that 
can serve not only as a weighting function for interpolation, but also as a smoothing kernel, as its effect is functionally similar to that of the PSF 
or any other function that can be applied to 
(see eq. [11]) before or during the sampling process leading to 
. It therefore follows that the weighting functions also possesses spectral filtering qualities, as the base spectrum 
is multiplied by its Fourier transform 
(see eq. [23]). One can, in fact, take advantage of this property in some cases. For example, the extraction of a signal from noise can be optimized by matching the spectral shape of the Fourier transform of the weighting function (more appropriately named the “filter” in this case) to that of the signal itself (if such information is available a priori). This is a result commonly established through the so‐called matched filter theorem (Haykin 1983). It is to be noted, however, that optimization of the S/N through filtering is not our goal. As is made evident in § 3.2, besides its fundamental role in the interpolation process, we are also concerned with determining the effects of the weighting function on the spatial resolution of a map. In general, the spatial extent of the smoothing kernel will always be significantly smaller than that of the optimized matched filter.
We now investigate the case of a map resulting from a resampling process where we seek to increase the density of samples. For example, a map with half the sampling periods as the original (i.e., of periods 
and 
) will consist of the combination of four different resampled maps [
, 
, 
, and 
] that all share the same sampling periods as the original map, but translated relative to each other. More precisely, we define 
In other words, 
is resampled at the same positions as the original map, while 
, 
, and 
are relatively shifted by 
, 
, and 
respectively. Calculating and summing the corresponding Fourier transforms (using eq. [23]), we find for the combined map that 
where 
is the normalization constant associated with 
. Correspondingly, we further define 
for 
, 3, 4, and we rewrite equation (29) as 
As is discussed presently in § 3.2, the magnitude of the coefficient 
depends on the width of the weighting function 
; the wider the function, the smaller the coefficient, and vice versa. As we see below, there are good reasons to limit the width of the weighting function, but if we assume for the moment that 
is such that 
, then 
and 
It is instructive to compare this result with the corresponding equation for an ESG 
, similar to that of equation (10), but with half the sampling interval in each direction: 
Again we see that apart from a multiplication factor, the (approximate) interpolated grid 
differs from 
by the presence of the convolution by 
. Although equations (31) and (32) are approximations and the aforementioned correspondence between 
and 
may fail for a given weighting function (i.e., 
in general in eq. [30]), this simplification is often reasonable (see § 3.2). At any rate, one can more closely approach the idealization of equation (32) by broadening the width of the weighting function 
(with a corresponding loss in spatial resolution, however).
3.2. Selection of the Weighting Function
There is a fair amount of subjectiveness in choosing the specific form and characteristics of a weighting function. We choose the following two criteria as guidelines for achieving this:
| 1. | The function must be sufficiently broad so that its amplitude is large enough at the interpolated positions, while not being too broad to significantly degrade the resolution of the map. | ||||
| 2. | Its spectral extent must be such that it filters out spatial frequencies for which | ||||
. Furthermore, we approximate the SHARC‐II PSF with the following Gaussian profile: 
The PSF size is usually defined by its full width at half‐maximum (FWHM), which is approximately 9
for SHARC‐II at 350 μm (Dowell et al. 2003). This gives 
; we will use standard deviations to specify widths of Gaussian PSFs. With this definition, the one‐sided bandwidth associated with SHARC‐II at 350 μm is 
Since 
, then 
and the Nyquist sampling criterion is met, as previously assumed. If we were to choose the weighting function 
to also be Gaussian and of width ϖ, then to satisfy criterion 2 above, we must have 
Taking the lower limit for the size of the kernel, we can evaluate the new resolution of the map with 
which corresponds to an equivalent PSF width of 9.6
for SHARC‐II and a loss of approximately 7% in spatial resolution. Finally, the relative amplitude of the weighting function at a distance of one‐half pixel away from the position of interpolation would be 
Although equations (35) and (36) satisfy criterion 1 above, and one could reasonably choose the corresponding weighting function to interpolate a map, one should nonetheless verify that the coefficients 
resulting from the interpolation process (see § 3.1) are sufficiently small when seeking to increase the density of samples in the final grid. Doing so will ensure that the approximation leading to equation (32) is adequate; for example. Figure 1 shows a map (top) of the normalization function 
for a SHARP ESG with the lower limit weighting function considered above (
). The top‐most curve in the lower part of the figure is a cut through a row or column of pixels for the normalization map. The bottom two curves are similar cuts for weighting functions of 
and 
, respectively. It is clear from these curves that the relative amplitude of the 
coefficients, which can be asserted from the level of ripple on the curves, exhibits a strong dependency on the width of the weighting function. For example, the two larger weighting functions in the lower part of Figure 1 exhibit amplitude variations of 11% and 1% for a loss in spatial resolution of 13% and 22%, respectively. This behavior is traced to the fact that a larger weighting function will more completely cover the space located between neighboring sampling positions; hence, the existence of a smoother normalization function 
. This will be better visualized for the general case with the graph shown in Figure 2, which plots trends in normalization function (solid line) and spatial resolution degradation (dashed lines) with smoothing kernel size (see the corresponding figure legend).
Fig. 1.— Map (top) of the normalization function 
for a SHARP ESG with a weighting function with 
(
for SHARP). The top‐most curve in the lower part of the figure is a cut through a row or column of pixels for the normalization map. The bottom two curves are similar cuts for weighting functions of 
and 
, respectively.
Fig. 2.— Trends in normalization function (solid line) and spatial resolution degradation (dashed lines) with smoothing kernel size. The kernel sizes on the abscissa are given as the Gaussian width (ϖ; bottom axis) and FWHM (= 
; top axis), both in units of the array pixel separation. The amplitude of the normalization function’s spatial variation is given as a percentage of the function’s average value (see Fig. 1). The corresponding loss in spatial resolution is described by eq. (35) and the following text. This is plotted here for different beam sizes and is indicated in units of pixel separation. For example, the case of SHARP, with 
, closely corresponds to the second dashed curve (from the top).
4. THE REGULARLY SAMPLED GRID
We define a regularly sampled grid (RSG) as being a generalization of the ESG discussed in the previous section. That is, a RSG has a well‐defined periodicity (just as the ESG), but it is a grid for which the pattern of Dirac distributions is more complex. While along a coordinate axis of the ESG there is only one Dirac distribution for a given period, a RSG may have many Dirac distributions (not necessarily evenly spaced) over the same interval. This difference is illustrated in cases a and b of Figure 3 for one‐dimensional versions of an ESG and a RSG, respectively. Practically speaking, a RSG would be encountered anytime that maps with similar characteristics, but that are translated relative to one another, are combined together to form a unique, final map. An example of this would be astronomical images of a given object at different pointing positions.
Fig. 3.— One‐dimensional examples of (a) an evenly sampled grid (ESG), (b) a regularly sampled grid (RSG), and (c)–(d) two irregularly sampled grids (ISGs). Interpolations at the positions of the two vertical broken lines would require weighting functions that have a common normalization factor for the ESG and RSG, but different normalization factors for the ISGs. Dirac distributions are shown as vertical arrows.
It should be apparent from this discussion and Figure 3b that a RSG, which we again denote by 
, can be simply expressed as a combination of a set of relatively displaced ESGs with 
where 
is defined in equation (11) and 
is the relative displacement associated with the pth of the 
ESGs that make up the RSG. It is straightforward to calculate the Fourier transform of equation (37) to get 
with 
as defined in equation (6).
The important aspect to emphasize for the interpolation of a RSG is that as was the case for an ESG, the normalization factor 
in equation (18) is common to all interpolated points similarly located within a one‐period segment anywhere on the grid. This is illustrated with the vertical broken lines in Figure 3. The existence of such a common normalization factor could be effectively adopted as the definition for a RSG.
If we interpolate our RSG using a weighting function 
that satisfies criterion 2 above, then the high spatial frequency components of the spectrum will be filtered out. Therefore, starting from equation (39), and using steps similar to those that led from equation (14) to equations (22) and (23), the resulting interpolated grid becomes 
and 
Once again c and 
respectively denote the common normalization factor and the displacement of the new interpolated grid in relation to the original grid (see eq. [21]). Obviously, the same comments apply for the map resulting from the resampling of a RSG here as for an ESG in § 3. That is, the most notable effect of the interpolation/resampling process is the presence of the convolution by the weighting function in equation (41).
5. THE IRREGULARLY SAMPLED GRID
An irregularly sampled grid (ISG) can manifest itself in different ways. For example, Figure 3c shows a case in which the distribution of Dirac functions within a given base period (of length l in the figure) is not the same from one interval to the next. Another possibility is shown in Figure 3d, where no Dirac distributions are present for some intervals. This can be likened to situations where pixels are missing from an array detector (see below).
The problem in the analysis of an ISG is twofold. First, there is no simple way of expressing the Fourier transform of irregularly spaced Dirac distributions such that the spectrum will show a repeating pattern of some frequency, as is the case for an ESG or a RSG. Moreover, the lack of regularity in the positions of the Dirac distributions implies that there does not exist a common normalization factor when performing interpolations to create an ESG from an ISG (see below). Nevertheless, it is still possible to analyze some specific types of ISGs. We deal with two possible cases in what follows.
5.1. The Combination of Relatively Translated and Rotated ESGs
It often happens that an astronomical source will be observed at different times, when it is at different locations and orientations on the celestial sphere. Invariably, we seek to combine the resulting images to form a final map of the object. If the array detector (which we assume to be perfect and therefore able to generate ESGs of data) used to record the images is part of an instrument that is unable to precisely track the apparent rotation of the source on the sky, then the different images of the source will be sampled with ESGs that will be rotated and possibly translated relative to each other. A simple example is shown in Figure 4, where two ESGs are combined: one rotated by 10° with respect to the other. These grids are not relatively translated (see the legend).
Fig. 4.— Combination of two relatively rotated ESGs. Every small filled circle corresponds to a Dirac distribution, and the position of the small open circle is the origin of the maps and of the rotation for one of the two grids. Its relative rotation is 10° with respect to the coordinate axes (also shown). Neither ESG is translated. The four large open circles correspond to the footprint of a predetermined weighting function.
Perhaps the most important aspect of Figure 4 is the fact that the combination of the two ESGs produces a grid that has an irregular pattern of Dirac distributions, as can be asserted by the coverage of a predetermined weighting function (shown by large open circles). Clearly, any weighting function 
is likely to cover a different number of samples at different positions on the map. The main consequence resulting from this fact will be the absence of a common normalizing factor at the different locations where interpolations are performed (note that the normalization function in eq. [19] will not be periodic). We can therefore expect that interpolated maps originating from ISGs will be more complex than those resulting from ESGs and RSGs.
Another point to consider is the possible relative rotation between the different maps to be combined. Because of this, we will do well to use the fact that the Fourier transform of a rotated map is the rotated Fourier transform of the original (i.e., without rotation) map. That is, if we have the Fourier pair 
then it is also true that 
(see the Appendix), where 
stands for the rotation operation (i.e., matrix). Because of this property of the Fourier transform, we can express an ISG 
composed of 
rotated and translated ESGs and its Fourier transform as 
where 
and 
are, respectively, the rotation matrix and the translation vector corresponding to grid p. Just as for the RSG, our goal is to generate an ESG 
of periods 
and 
(along the x‐ and y‐axes, respectively) from 
. Although, as was previously pointed out, we cannot express the interpolation process with a simple convolution with a weighting function 
, we can still use the general expression given in equation (20). That is, with 
denoting the origin of the new interpolated grid (see eq. [21]), we have 
Calculating the Fourier transform of equation (44), we get 
which, using the assumption that 
is such that it filters out the higher frequency components of the spectrum, can be approximated to 
Correspondingly, we can approximate equation (44) to 
Equation (45) shows best the effect of the lack of a common normalization factor on the interpolation process. Since the Fourier transform 
of the normalization function is convolved with the weighted (and low‐pass–filtered) spectrum 
, the resulting spectrum of the interpolated ESG 
is broadened by 
. It is interesting to note that 
and 
have opposite effects on the signal. That is, the weighting function restricts the extent of the spectrum, while the normalization function extends it.
Given such an ISG, it should in principle be possible to evaluate 
and quantify its effect. In particular, one should ensure that the two previous criteria (see § 3.2) used to select the weighting function are met. Optimally, 
will exhibit slow variations and sufficiently low amplitude that the spectral broadening will be minimal. To make this clearer, we show in Figure 5 an example consisting of a combination of two relatively rotated ESGs, similar to those of Figure 4 (i.e., one rotated by 10° with respect to the other, and no relative translation between the two). The map at the top of the figure is for the normalization function 
of the resulting SHARP ISG using a weighting function with 
(
for SHARP). The top‐most curve in the lower part of the figure is an arbitrary cut through a row of pixels for this normalization map. The bottom two curves are similar cuts for weighting functions of 
and 
, respectively. The black dots highlight the values taken by 
for a resampling onto an ESG at the original sampling rate. It is clear from the top two cuts that the normalization factor is not constant in general. One can also assert from this that the spectrum due to a broad source (in relation to the size of the map) would be significantly more broadened by the weighting function that produced the top curve (i.e., with 
) than by the other two.
Fig. 5.— Combination of two relatively rotated ESGs, similar to those of Fig. 4. We show a map (top) of the normalization function 
for the resulting SHARP ISG using a weighting function with 
(
for SHARP). The top‐most curve in the bottom part of the figure is an arbitrary cut through a row of pixels for this normalization map. The bottom two curves are similar cuts for weighting functions of 
and 
, respectively. The filled circles highlight the values taken by 
for a resampling onto an ESG at the original sampling rate. It is clear from the top two cuts that the normalization factor is not constant in general.
5.2. The Effects of Missing Samples
It is a common, if unfortunate, fact that detector arrays used in astronomy will often contain pixels that are either performing significantly below specifications or are completely unusable. Astronomers usually work around the difficulties occasioned by these so‐called missing pixels by dithering the array during observations, thus ensuring complete mapping of the source under study. However, it would be instructive to analyze and quantify the impact that missing pixels would have on the representation of astronomical signals without such corrective techniques.
Although we take into account the fact that the map obtained from the array is composed of a finite number of samples, our approach consists of first temporarily lending it an infinite character and then removing it. More precisely, although the size of the (rectangular) detector array considered in this section is 
, we first assume that the two‐dimensional pattern of 
pixels (including the missing pixels) repeats infinitely in all directions. The underlying assumption is that the finiteness of the map will be restored in the end by windowing with the appropriate aperture function. Using this approach, we express the sampled signal (before windowing) as 
where 
is the position of a pixel on the array and (take note of the periods) 
That is, we first account for the pixels of the finite array through the summation 
, and then associate a Dirac train of periods 
to each pixel. These Dirac trains are relatively translated in space and are accounted for by the summations on s and t in equation (47). The Fourier transform of the combined Dirac trains is 
with 
Note that the minimum separation between two Dirac distributions in frequency space is 
, where Nl is the greater of 
and 
. We write the right‐hand side of equation (49), which we denote by 
, as follows: 
with 
The effect of missing pixels can now easily be quantified, at least in principle, by omitting them from the 
summation in equation (52). As an example, consider the case in which the 
pixels in the “top right” corner of the array are missing. For this particular example, we will do well to make the following substitutions: 
where the indices i and k stand for the columns and rows of the detector, respectively. We can then transform equation (52) to 
When 
(i.e., when all the pixels are accounted for), equation (53) simplifies to 
where 
and 
are some integer numbers. That is to say, equation (51) then becomes 
which is, as it should be, the same result as was obtained earlier with equation (4) for the ESG.
Although it is necessary to plot 
to assess the effect of an arbitrary distribution of missing pixels, it should now be clear from equation (53) that its amplitude is in general nonzero for all values of m and n. In other words, by removing even only one pixel, we went from a case in which we only had Dirac functions at frequency intervals of 
and 
to a situation where they are separated by intervals of only 
and 
. It is important, however, to quantify the relative magnitude of these Dirac distributions. For example, returning to equation (53) pertaining to our case of the 
missing “top right” pixels, we find that 
This last relation yields the perhaps intuitive result that when only a small number of pixels are missing, the amount of contamination determined by the ratio expressed in equation (54) is approximately equal to the ratio of the number of missing pixels to the number of good pixels. However, we should resist the temptation to generalize this result, since different distributions of 
missing pixels would give different levels of contamination; especially if they are not concentrated in one part of the array, as is the case here. An example is shown in Figure 6, where the function 
is plotted for three different cases. Starting with the nominal SHARP 
array, 
is shown for 
when (1) no pixels (black curve and dots), (2) 16 randomly positioned pixels (red curve and dots), and (3) the 
“top right” corner pixels (blue curve and dots) are missing. Contrary to the case of an ESG (corresponding to the black dots) where 
when 
, an ISG (red and blue dots) will in general have 
for all m and n. It should be clear that 
acts as a mask that will or will not allow the appearance of Dirac distributions that are more closely spaced in frequency, depending on whether or not there are missing pixels in the array. This serves to emphasize the fact that missing pixels will bring some spectral contamination (i.e., aliasing) in the sampled signal.
Fig. 6.— Function 
for 
when no pixels (black), 16 randomly positioned pixels (red), and the 
“top right” corner pixels (blue) are missing. Contrary to the case of an ESG (corresponding to the black dots) where 
when 
, an ISG (red and blue dots) will in general have 
for all m and n. The amplitude of 
for the relevant (i.e., integer) values of m are shown by the colored dots. The curves, which include computations at intermediate values of m, are only shown to emphasize the fact that 
can be interpreted as a mask function (see text).
This becomes more evident if we calculate the spectrum of the measured signal from equations (47) and (51): 
We see that the different replicas of 
are spaced in frequency according to equation (50) (compare this with the case of the ESG in eq. [6], where the spacing between replicas is N times greater). Contrary to the case of an ESG, a convolution with the usual weighting function 
will not restore a low‐pass–filtered version of the 
map. To make this clear, we first define a new function 
such that 
Then proceeding with the usual interpolation defined in equation (18), we find that 
and 
where 
is defined in equation (21) and denotes, once again, the position of the origin of the interpolated grid. These equations show that maps resulting from the interpolation of ISGs containing missing pixels suffer from both spectral aliasing [from the presence of 
in lieu of 
in eq. (57)] and broadening [because of the presence of 
].
We once again stress the realization that missing pixels will bring some amount of aliasing that will be impossible to remove with a reasonably sized weighting function. The concept of Nyquist sampling can even lose much of its meaning and usefulness in a situation where too many pixels are missing, or when the level of contamination due to spectral aliasing is comparable to the noise level present in the map . This fact strongly underlines the necessity of performing adequate dithers or other scanning strategies when observing with an imperfect detector array.
We complete the analysis by performing the windowing mentioned earlier, which transforms equation (55) to 
For reasonably large detector arrays, we do not expect that the presence of the sinc functions will be of any significance, due to their spectral narrowness relative to the extent of 
[or 
]. Because of this, our periodic depiction of the array, on which our analysis rests, is justified.
6. SUMMARY
In this paper, we addressed the question of astronomical image processing from data obtained with array detectors. We defined and analyzed the cases of evenly (ESG), regularly (RSG), and irregularly (ISG) sampled grids for idealized and realistic detectors. We focused on the effect of interpolation on the maps, while using a Gaussian kernel to accomplish this task. In all cases (i.e., ESG, RSG, and ISG), we have applied the method of weighted averages (eq. [18]) to produce a map interpolated on a finely spaced grid.
We defined an ESG as a map where the signal to be analyzed is digitized with a simple, two‐dimensional train of Dirac distributions evenly separated with well‐defined spacings (see eq. [3]). Moreover, since the ESG is the simplest way to represent and analyze a set of sampled data, we always sought to transform a non‐evenly sampled grid (i.e., a RSG or an ISG) to an ESG through the process of interpolation. While studying the ESG, we found that the interpolation process invariably leads to a loss in spatial resolution, and that this loss grows with increasing width of the smoothing kernel (this result is true in general; i.e., when considering RSGs and ISGs). When an ESG is resampled at the same rate as the original map, the interpolation process can usually be adequately taken into account by replacing the original signal by its convolution with the weighting function (see eq. [24]). However, the same is not true in general when the final ESG is resampled at a rate different than that of the original map (see eq. [30]). This is due to the fact that the normalization function that is intrinsic (and necessary) to the interpolation process is not constant but is a function of position (although it is periodic, with periods corresponding to the sampling rates). The aforementioned replacement of the original map in the interpolation process by its convolution with the weighting function will only be adequate in such cases when the latter is sufficiently broad relative to the spacing between samples (see Figs. 1 and 2).
We defined a RSG as a generalization of an ESG such that it consists of a combination of a number of relatively translated ESGs of similar sampling rates. All the results obtained for the ESG can be generalized to the RSG.
We analyzed two different types of ISGs: the combination of relatively translated and rotated ESGs, and ESG‐like maps with missing samples (e.g., data grids made with detectors exhibiting dead pixels). In the first case, the interpolation process cannot be simply represented by a convolution of the original signal with the weighting function, as this operation must be subsequently multiplied by the normalization function (see eq. [46]). Because of the irregular nature of the new map, the normalization function is not periodic in general and will add structure to the spectrum (i.e., the Fourier transform) of the map. More precisely, the spectrum of the source (filtered by the spectral profile of the weighting function) will be broadened through its convolution with the Fourier transform of the normalization function. This effect is a function of the size of the weighting function, as the spatial variation of the normalization function grows larger for smaller kernel widths. Although these results also apply to maps exhibiting missing samples, we found that these further suffer from spectral aliasing that may reduce or negate the usefulness of the Nyquist sampling criterion in extreme cases. This fact strongly underlines the necessity of performing adequate dithers or other scanning strategies when observing with an imperfect detector array.
M. H.’s research is funded through the NSERC Discovery Grant, Canada Research Chair, Canada Foundation for Innovation, Ontario Innovation Trust, and Western’s Academic Development Fund programs. J. E. V. acknowledges support from NSF grants AST 05‐40882 to the California Institute of Technology and AST 05‐05124 to the University of Chicago. SHARC II is also funded through the NSF grant AST 05‐40882 to the California Institute of Technology. SHARP is funded through the NSF grants AST 02‐43156 and AST 05‐05230 to Northwestern University.
APPENDIX A
In this Appendix, we provide a few simple derivations to justify some of the results used in the text.
A1. FOURIER TRANSFORM OF A DIRAC TRAIN
The Fourier transform of a Dirac train can easily be determined by first calculating the associated Fourier series. In the one‐dimensional case, we have 
since the Fourier series for a periodic function 
of period l is defined as 
where the Fourier coefficient 
is 
Before calculating the Fourier transform of equation (A1), we consider the so‐called duality property of the Fourier pair in equations (1) and (2). More precisely, we mean that if for a function 
we have 
then it must also be true that for 
we have 
as can be readily verified by inspection of equations (1) and (2). It follows from this that since 
then 
Using the one‐dimensional version of equation (A2) to calculate the Fourier transform of equation (A1), we find that 
The two‐dimensional generalization of this result is straightforward and leads to equations (3) and (4).
A2. THE NYQUIST SAMPLING CRITERION
The “Nyquist sampling criterion” can be understood with equation (A3) and the product/convolution property of the Fourier transform. This property states that the following Fourier pair is valid 
for two functions 
and 
, as can easily be verified from the definition of the Fourier transform. Therefore, for the sampling of a function 
, we have 
which implies that the base spectrum 
is repeated in frequency space at an interval equal to the sampling period of 
. It is apparent that in order to avoid any cross‐contamination between the different spectral replicas of 
, the following relation must be enforced: 
Relation (A4) is the one‐dimensional mathematical equivalent of the Nyquist sampling criterion, which states that in order to recover the base spectrum 
from its sampled version (through spectral filtering), the sampling frequency must be at least twice as large as the frequency extent of 
.
A3. NORMALIZATION FACTOR
We know from our analysis that an interpolated map 
resulting from a previously sampled data set 
is given by 
which we transform slightly to 
As usual, 
and 
are the normalization and weighting functions, respectively, and the vectors 
and 
are given by equations (5) and (21). When the original map is an ESG, the normalization function is periodic and can therefore be expanded with a two‐dimensional Fourier series 
where 
is the corresponding Fourier coefficient and 
is given by equation (6).
From equations (A2) and (A6), we can write the Fourier transform of equation (A5) as 
but since 
with 
and 
, then 
However, we can use equation (A6) one more time to transform this last relation to 
and 
Evidently, 
is constant for a given interpolated map, but will vary as a function of the displacement of the new sampling grid (defined by 
) relative to the original one. Equation (A7) leads to (and justifies) equation (24), provided that we set 
.
A4. FOURIER TRANSFORM OF A ROTATED MAP
Finally, we prove the result used in § 5 that the Fourier transform of a rotated map is the rotated Fourier transform of the unrotated map. To do so, we subject a two‐dimensional map 
to a rotation 
and calculate the Fourier transform of the transformed map 
with 
We now make the change of variable 
to get 
where the last transformation was made possible by the fact that the inverse of a rotation matrix equals its transpose. We therefore find from equation (A8) that if 
then 
REFERENCES
- Dowell, C. D., et al. 2003, Proc. SPIE 4855, 73
- Granrath, D., & Lersch, J. 1998, J. Opt. Soc. Am. A, 15, 791
- Haykin, S. S. 1983, Communication Systems (2nd ed.; New York: Wiley)
- Lombardi, M., & Schneider, P. 2001, A&A, 373, 359
- Novak, G., et al. 2004, Proc. SPIE 5498, 278
-
1 In this paper, we use a lowercase letter and the corresponding capital letter for a function and its Fourier transform, respectively.
-
2 It should be noted that the normalization function will be constant for an infinite grid when
for 
or 
(e.g., sinc weighting functions of corresponding widths in normal space). This condition must be strictly enforced in order to obtain a constant normalization factor while satisfying the Nyquist sampling criterion.