2.9 Moiré pattern in stacking and its correction ¶

After warping some images with the default mode of Warp (see Align pixels with WCS considering distortions) you may notice that the background noise is no longer flat. Some regions will be smoother and some will be sharper; depending on the orientation and distortion of the input/output pixel grids. This is due to the Moiré pattern, which is especially noticeable/significant when two slightly different grids are super-imposed.

With the commands below, we’ll download a single exposure image from the J-PLUS survey and run Warp (on a \(8\times8\) arcmin\(^2\) region to speed it up the demos here). Finally, we’ll open the image to visually see the artificial Moiré pattern on the warped image.

## Download the image (73.7 MB containing an 9216x9232 pixel image)
$ jplusdr2=http://archive.cefca.es/catalogues/vo/siap/jplus-dr2/reduced
$ wget $jplusdr2/get_fits?id=771463 -Ojplus-exp1.fits.fz

## Align a small part of it with the sky coordinates.
$ astwarp jplus-exp1.fits.fz --center=107.62920,39.72472 \
          --width=8/60 -ojplus-e1.fits

## Open the aligned region with DS9
$ astscript-fits-view jplus-e1.fits

In the opened DS9 window, you can see the Moiré pattern as wave-like patterns in the noise: some parts of the noise are more smooth and some parts are more sharp. Right in the center of the image is a blob of sharp noise. Warp has the --checkmaxfrac option for direct inspection of the Moiré pattern (described with the other options in Align pixels with WCS considering distortions). When run with this option, an extra HDU (called MAX-FRAC) will be added to the output. The image in this HDU has the same size as the output. However, each output pixel will contain the largest (maximum) fraction of area that it covered over the input pixel grid. So if an output pixel has a value of 0.9, this shows that it covered \(90\%\) of an input pixel. Let’s run Warp with --checkmaxfrac and see the output (after DS9 opens, in the “Cube” window, flip between the first and second HDUs):

$ astwarp jplus-exp1.fits.fz --center=107.62920,39.72472 \
          --width=8/60 -ojplus-e1.fits --checkmaxfrac

$ astscript-fits-view jplus-e1.fits

By comparing the first and second HDUs/extensions, you will clearly see that the regions with a sharp noise pattern fall exactly on parts of the MAX-FRAC extension with values larger than 0.5. In other words, output pixels where one input pixel contributed more than half of the its value. As this fraction increases, the sharpness also increases because a single input pixel’s value dominates the value of the output pixel. On the other hand, when this value is small, we see that many input pixels contribute to that output pixel. Since many input pixels contribute to an output pixel, it acts like a convolution, hence that output pixel becomes smoother (see Spatial domain convolution). Let’s have a look at the distribution of the MAX-FRAC pixel values:

$ aststatistics jplus-e1.fits -hMAX-FRAC
Statistics (GNU Astronomy Utilities) 0.23
-------
Input: jplus-e1.fits (hdu: MAX-FRAC)
-------
  Number of elements:                      744769
  Minimum:                                 0.250213461
  Maximum:                                 0.9987495374
  Mode:                                    0.5034223567
  Mode quantile:                           0.3773819498
  Median:                                  0.5520805544
  Mean:                                    0.5693956458
  Standard deviation:                      0.1554693738
-------
Histogram:
 |                      ***
 |                   **********
 |                 *****************
 |              ************************
 |           *******************************
 |         **************************************
 |       *********************************************
 |     ****************************************************
 |   ***********************************************************
 | ******************************************************************
 |**********************************************************************
 |----------------------------------------------------------------------

The smallest value is 0.25 (=1/4), showing that 4 input pixels contributed to the output pixels value. While the maximum is almost 1.0, showing that a single input pixel defined the output pixel value. You can also see that the most probable value (the mode) is 0.5, and that the distribution is positively skewed.

This is a well-known problem in astronomical imaging and professional photography. If you only have a single image (that is already taken!), you can undersample the input: set the angular size of the output pixels to be larger than the input. This will decrease the resolution of your image, but will ensure that pixel-mixing will always happen. In the example below we are setting the output pixel scale (which is known as CDELT in the FITS standard) to \(1/0.5=2\) of the input’s. In other words each output pixel edge will cover double the input pixel’s edge on the sky, and the output’s number of pixels in each dimension will be half of the previous output.

$ cdelt=$(astfits jplus-exp1.fits.fz --pixelscale -q \
                  | awk '{print $1}')
$ astwarp jplus-exp1.fits.fz --center=107.62920,39.72472 \
          --width=8/60 -ojplus-e1.fits --cdelt=$cdelt/0.5 \
          --checkmaxfrac

In the first extension, you can hardly see any Moiré pattern in the noise. When you go to the next (MAX-FRAC) extension, you will see that almost all the pixels have a value of 1. Of course, decreasing the resolution by half is a little too drastic. Depending on your image, you may be able to reach a sufficiently good result without such a drastic degrading of the input image. For example, if you want an output pixel scale that is just 1.5 times larger than the input, you can divide the original coordinate-delta (or “cdelt”) by \(1/1.5=0.6666\) and try again. In the MAX-FRAC extension, you will see that the range of pixel values is now between 0.56 to 1.0 (recall that originally, this was between 0.25 and 1.0). This shows that the pixels are more similarly mixed and in fact, when you look at the actual warped image, you can hardly distinguish any Moiré pattern in the noise.

However, deep astronomical data are usually built by several exposures (images), not a single one. Each image is also taken by (slightly) shifting the telescope compared to the previous exposure. This shift is known as “dithering” or a “pointing pattern”, see Pointing pattern design. We do this for many reasons (for example tracking errors in the telescope, high background values, removing the effect of bad pixels or those affected by cosmic rays, robust flat pattern measurement, etc.⁷⁸). One of those “etc.” reasons is to correct the Moiré pattern in the final coadded deep image.

The Moiré pattern is fixed to the grid of the image, slightly shifting the telescope will result in the pattern appearing in different parts of the sky. Therefore when we later stack, or coadd, the separate exposures into a deep image, the Moiré pattern will be decreased there. However, dithering has possible drawbacks based on the scientific goal. For example when observing time-variable phenomena where cutting the exposures to several shorter ones is not feasible. If this is not the case for you (for example in galaxy evolution), continue with the rest of this section.

Because we have multiple exposures that are slightly (sub-pixel) shifted, we can also increase the spatial resolution of the output. For example, let’s set the output coordinate-delta (--cdelt, or pixel scale) to be 1/2 of the input. In other words, the number of pixels in each dimension of the output is double the first Warp command of this section:

$ astwarp jplus-exp1.fits.fz --center=107.62920,39.72472 \
          --width=8/60 -ojplus-e1.fits --cdelt=$cdelt/2 \
          --checkmaxfrac

$ aststatistics jplus-e1.fits -hMAX-FRAC --minimum --maximum
6.26360438764095e-02 2.50680270139128e-01

$ astscript-fits-view jplus-e1.fits

From the last command, you see that like the previous change in --cdelt, the range of MAX-FRAC has decreased. However, when you look at the warped image and the MAX-FRAC image with the last command, you still visually see the Moiré pattern in the noise (although it has significantly decreased compared to the original resolution). It is still present because 2 is an exact multiple of 1. Let’s try increasing the resolution (oversampling) by a factor of 1.25 (which isn’t an exact multiple of 1):

$ astwarp jplus-exp1.fits.fz --center=107.62920,39.72472 \
          --width=8/60 -ojplus-e1.fits --cdelt=$cdelt/1.25 \
          --checkmaxfrac
$ astscript-fits-view jplus-e1.fits

You don’t see any Moiré pattern in the noise any more, but when you look at the MAX-FRAC extension, you see it is very different from the ones you had seen before. In the previous MAX-FRAC image, you could see large blobs of similar values. But here, you see that the variation is almost on a pixel scale, and the difference between one pixel to the next is not significant. This is why you don’t see any Moiré pattern in the warped image.

In J-PLUS, each part of the sky was observed with a three-point pointing pattern (very small shifts in each pointing). Let’s download the other two exposures and warp the same region of the sky to the same pixel grid (using the --gridfile feature). Then, let’s open all three warped images in one DS9 instance:

$ wget $jplusdr2/get_fits?id=771465 -Ojplus-exp2.fits.fz
$ wget $jplusdr2/get_fits?id=771467 -Ojplus-exp3.fits.fz

$ astwarp jplus-exp2.fits.fz --gridfile jplus-e1.fits \
          -o jplus-e2.fits --checkmaxfrac
$ astwarp jplus-exp3.fits.fz --gridfile jplus-e1.fits \
          -o jplus-e3.fits --checkmaxfrac

$ astscript-fits-view jplus-e*.fits

In the three warped images, you don’t see any Moiré pattern, so far so good... now, take the following steps:

1. In the small “Cube” window, click the “Next” button so you see the MAX-FRAC extension/HDU.
2. Click on the “Frame” button (in the top row of buttons just on top of the image), and select the “Single” button in the bottom row.
3. Open the “Zoom” menu (not button), and select “Zoom 16”.
4. Press the TAB key to flip through each exposure.
5. Focus your eyes on the pixels with the largest value (white colored pixels), while pressing TAB to flip between the exposures. You will see that in each exposure they cover different pixels (nicely getting averaged out after stacking).

The exercise above shows that the Moiré pattern (that had already decreased significantly) will be further decreased after we stack the images. So let’s stack these three images with the commands below. First, we need to remove the sky-level from each image using NoiseChisel, then we’ll stack the INPUT-NO-SKY extensions using filled MAD-clipping (to reject outliers, and especially diffuse outliers, robustly, see Clipping outliers).

$ for i in $(seq 3); do \
   astnoisechisel jplus-e$i.fits -ojplus-nc$i.fits; \
  done

$ astarithmetic jplus-nc*.fits 3 5 0.2 sigclip-mean \
                -gINPUT-NO-SKY -ojplus-stack.fits

$ astscript-fits-view jplus-nc*.fits jplus-stack.fits

After opening the individual exposures and the final stack with the last command, take the following steps to see the comparisons properly:

1. Click on the stack image so it is selected.
2. Go to the “Frame” menu, then the “Lock” item, then activate “Scale and Limits”.
3. Scroll your mouse or touchpad to zoom into the image.

You clearly see that the stacked image is deeper and that there is no Moiré pattern, while you have slightly improved the spatial resolution of the output compared to the input. In case you want the stack to have the original pixel resolution, you just need one more warp:

$ astwarp jplus-stack.fits --cdelt=$cdelt -ojplus-stack-origres.fits

For optimal results, the oversampling should be determined by the dithering pattern of the observation: For example if you only have two dither points, you want the pixels with maximum value in the MAX-FRAC image of one exposure to fall on those with a minimum value in the other exposure. Ideally, many more dither points should be chosen when you are planning your observation (not just for the Moiré pattern, but also for all the other reasons mentioned above). Based on the dithering pattern, you want to select the increased resolution such that the maximum MAX-FRAC values fall on every different pixel of the output grid for each exposure. Note that this discussion is on small shifts between pointings (dithers), not large ones like offsets); see Pointing pattern design.