A Stacking Case Study
First we are going to load 8 frames of data that is included with the
ProPane package:
image_list = Rfits_make_list(dirlist = system.file('extdata/stack/', package="ProPane"),
extlist = 2) #extlist=2 because these are compressed images
In general the above is how you will probably load data for stacking.
You can also fully load the images if the resultant list is small enough
(easily fits within memory). That will be a bit faster in cases where it
is possible.
Let’s check the loaded image pointers. They are overlapping regions
of sky:
for(i in 1:8){
plot(image_list[[i]])
}
Now we will create a target WCS for stacking (if not provided then
the WCS of the first image is used instead):
keyvalues_out = Rwcs_setkeyvalues(
CRVAL1 = 36.8962, #Roughly the middle of the images
CRVAL2 = -5.1906, #Roughly the middle of the images
pixscale = 0.3, #this is a bit finer than the input pixel scale of 0.337
NAXIS1 = 2200,
NAXIS2 = 2200
)
Now we make a simple stack:
stack = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out,
magzero_in = 30,
magzero_out = 23.9 #micro-jansky output
)
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Time taken: 9.731 seconds
plot(stack$image, qdiff=TRUE)
plot(stack$weight, magmap=FALSE)
We can see how deep most of the image is:
table(as.integer(stack$weight$imDat)) / prod(dim(stack$weight))
#>
#> 0 1 2 3 4 5
#> 0.083463223 0.018859091 0.007114669 0.008506818 0.012958264 0.008068595
#> 6 7 8
#> 0.013600620 0.015836777 0.831591942
So about 83% of the image has 8 frames contributing.
We can visually see the change in depth:
propaneWarp(image_list[[1]][,], keyvalues_out = keyvalues_out, plot=TRUE, qdiff=TRUE)
We can also check the pixel statistics of the stack:
par(mar=c(3.1,3.1,1.1,1.1))
maghist(stack$image$imDat[stack$weight$imDat==8], xlim=5)
#> Summary of used sample:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.1919930 -0.0262250 -0.0003756 0.0010021 0.0262025 0.1950983
#> Pop Std Dev: 0.041015
#> MAD: 0.038848
#> Half 16-84 Quan (1s): 0.039097
#> Half 02-98 Quan (2s): 0.082324
#> Using 3974588 out of 4024905 (98.75%) data points (9 < xlo & 50308 > xhi)
This implies the sky RMS in the deepest region is about 0.039 ADU
(which happens to be in micro-jansky). We can compare this to the
predicted:
1/sqrt(median(stack$inVar$imDat[stack$weight$imDat==8]))
#> [1] 0.04583811
This is a bit larger at about 0.046 ADU. The main source of the
difference is pixel covariance, which has the effect of artificially
reducing the apparent RMS.
A Coarser Case
Now we are going to look at using a coarser pixel scale on the target
WCS:
keyvalues_out_2 = Rwcs_setkeyvalues(
CRVAL1 = 36.8962, #Roughly the middle of the images
CRVAL2 = -5.1906, #Roughly the middle of the images
pixscale = 1, #this is coarser than the input pixel scale of 0.337
NAXIS1 = 800,
NAXIS2 = 800
)
First we will do a forward projection, where all source pixels are
caste onto the target WCS (which is like drizzling):
stack_2_forward = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out_2,
magzero_in = 30,
magzero_out = 23.9, #micro-jansky output
direction = 'forward'
)
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Time taken: 6.2 seconds
Then we can try back projection, where we use the final WCS to look
up input fluxes (more like Swarp):
stack_2_backward = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out_2,
magzero_in = 30,
magzero_out = 23.9, #micro-jansky output
direction = 'backward'
)
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Time taken: 7.256 seconds
And now let’s look at them:
plot(stack_2_forward$image, qdiff=TRUE)
plot(stack_2_backward$image, qdiff=TRUE)
Is there a different in apparent depth?
par(mar=c(3.1,3.1,1.1,1.1))
maghist(stack_2_forward$image$imDat, xlim=5)
#> Summary of used sample:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -1.331391 -0.182112 -0.010753 0.004158 0.168800 1.311032
#> Pop Std Dev: 0.29366
#> MAD: 0.25997
#> Half 16-84 Quan (1s): 0.2644
#> Half 02-98 Quan (2s): 0.60931
#> Using 393249 out of 640000 (61.45%) data points (258 < xlo & 7848 > xhi)
par(mar=c(3.1,3.1,1.1,1.1))
maghist(stack_2_backward$image$imDat, xlim=5)
#> Summary of used sample:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -2.253146 -0.301089 -0.002096 0.013412 0.303904 2.255121
#> Pop Std Dev: 0.48724
#> MAD: 0.44823
#> Half 16-84 Quan (1s): 0.45148
#> Half 02-98 Quan (2s): 0.98609
#> Using 395468 out of 640000 (61.79%) data points (211 < xlo & 4898 > xhi)
Indeed- the forward propagated data is deeper, which makes sense
because more pixels are caste into each target pixel. The backwards
stack is still deeper than a single frame since we have still stacked 8
images! The forward propagated data is however much more covariant:
cor(as.numeric(stack_2_forward$image$imDat[1:799,1:799]),
as.numeric(stack_2_forward$image$imDat[2:800,2:800]), use='complete.obs')
#> [1] 0.4103285
cor(as.numeric(stack_2_backward$image$imDat[1:799,1:799]),
as.numeric(stack_2_backward$image$imDat[2:800,2:800]), use='complete.obs')
#> [1] 0.2898906
The default is to project from the higher resolution grid, so in this
case without specifiying Rwcs_stack would carry out
forward propagation. Which makes more sense is really a choice for the
user.
A Quick Clip
Commonly when we have multiple input frames we want to detect outlier
hot/cold pixels and clip them out from the stacked image. In
Rwcs_Stack this is achieved with the doclip
option. To speed up the re-stacking of frames we can set
dump_frames to TRUE (meaning warped frames are saved to a
temporary directory and loaded back in as needed, saving re-projection
time when ignoring clipped pixels). To see information regarding
clipping we can also set keep_extreme_pix to TRUE.
stack_clip = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out,
magzero_in = 30,
magzero_out = 23.9, #micro-jansky output
keep_extreme_pix = TRUE,
keepcrop = FALSE,
doclip = TRUE,
dump_frames = TRUE,
return_all = TRUE
)
#> Frames being dumped to /var/folders/0j/jjlx9l6d6j75ffvr5l7mcg7w0000gn/T//RtmpNOvfYn
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Calculating Extreme Pixels 1 to 8 of 8
#> Clipping out extreme cold/hot pixels
#> Restacking without clipped cold/hot pixels
#> Stacking Images and InVar 1 to 8 of 8
#> Time taken: 19.94 seconds
The resulting stacking looks very similar, but we can see small
differences:
plot(stack$image - stack_clip$image)
#> [1] "Too many same valued pixels: turning off magmap scaling!"
We can check the actual clipped pixels directly:
plot(stack_clip$clip)
#> [1] "Too many same valued pixels: turning off magmap scaling!"
And we can zoom in to see the difference in more detail:
plot(stack$image[1000,500,box=500])
plot(stack_clip$image[1000,500,box=500])
In this case the clipped sources look a bit odd. We saved the input
warped frames so we can check out all 8:
for(i in 1:8){
plot(stack_clip$image_pre_stack[[i]][1000,500,box=500])
}
Frames 1 and 4 seem to have unusual sources (perhaps cosmic rays)
that are not present in the other frames. It is notable that our stack
has not been aggressive enough in removing these bad regions, so we can
tweak the clip_tol:
stack_clip2 = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out,
magzero_in = 30,
magzero_out = 23.9, #micro-jansky output
keep_extreme_pix = TRUE,
keepcrop = FALSE,
doclip = TRUE,
clip_tol = 10,
dump_frames = FALSE,
return_all = TRUE
)
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Calculating Extreme Pixels 1 to 8 of 8
#> Clipping out extreme cold/hot pixels
#> Restacking without clipped cold/hot pixels
#> Stacking Images and InVar 1 to 8 of 8
#> Time taken: 17.63 seconds
plot(stack_clip2$clip[1000,500,box=500])
#> [1] "Too many same valued pixels: turning off magmap scaling!"
plot(stack_clip2$image[1000,500,box=500])
That looks much better! The clip_tol is effectively the
sigma clipping to apply, so setting it to a lower number like 5-10 is
often sensible (the default is a deliberately very high 100).
Another route is to use the more conservative clip_tol but
specify that we dilate around the clipped pixels (to capture bleeding
around cosmic rays etc):
stack_clip3 = propaneStackWarpInVar(image_list = image_list,
inVar_list = 0.000468, #assume all have the same inVar
exp_list = 10, #all have the same exposure time of 10s
cores = 8,
keyvalues_out = keyvalues_out,
magzero_in = 30,
magzero_out = 23.9, #micro-jansky output
keep_extreme_pix = TRUE,
keepcrop = FALSE,
doclip = TRUE,
clip_tol = 50,
clip_dilate = 5,
return_all = TRUE #This will save the warped outputs
)
#> Stacking 8 images; Nbatch: 8; cores: 8
#> Projecting Images 1 to 8 of 8
#> Projecting Inverse Variance 1 to 8 of 8
#> Projecting Exposure Times 1 to 8 of 8
#> Stacking Images and InVar 1 to 8 of 8
#> Stacking Exposure Times 1 to 8 of 8
#> Calculating Extreme Pixels 1 to 8 of 8
#> Clipping out extreme cold/hot pixels
#> Restacking without clipped cold/hot pixels
#> Stacking Images and InVar 1 to 8 of 8
#> Time taken: 17.95 seconds
plot(stack_clip3$clip[1000,500,box=500])
#> [1] "Too many same valued pixels: turning off magmap scaling!"
plot(stack_clip3$image[1000,500,box=500])
The best route to clipping is very much data dependent, so some
experimentation from the user might be required. In general you will
want to set in the range 5-100 though.
Note if you are doing clipping on a lot of frames then it will not be
feasible to run with return_all = TRUE since this will keep all
the warped projections in memory. The better option in that case is to
run with dump_frames = TRUE, which will save the warped frames
to a temporary directory and load them as needed when applying the
per-frame clipping.
A Median Clip
You can also take the warped frames to compute a median stack. The
advantage of this can be robustness to very bad data, but the caveat is
that in the realm of well behaved data a median stack has an RMS larger
by a factor \(\sqrt{\pi/2} = 1.253\)
(no such thing as money for nothing). That might not sound like a lot,
but it basically means you would need to observe everything 50% longer
again to get the same signal-to-noise when doing median stack (which is
quite a big overhead).
library(imager)
#> Loading required package: magrittr
#>
#> Attaching package: 'imager'
#> The following object is masked from 'package:magrittr':
#>
#> add
#> The following objects are masked from 'package:stats':
#>
#> convolve, spectrum
#> The following object is masked from 'package:graphics':
#>
#> frame
#> The following object is masked from 'package:base':
#>
#> save.image
stack_med = propaneStackFlatFunc(stack_clip3$image_pre_stack, imager_func = parmed)
magimage(stack_med$image); legend('topleft', legend='Median Stack')
magimage(stack_clip3$image$imDat); legend('topleft', legend='inVar Stack')
magimage(stack_med$image - stack_clip3$image$imDat, qdiff=TRUE); legend('topleft', legend='Diff')
We can see there is a difference in depth when comparing the image
histograms:
par(mar=c(3.1,3.1,1.1,1.1))
maghist(stack_med$image, xlim=5)
#> Summary of used sample:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.2277246 -0.0302569 -0.0001985 0.0012074 0.0306361 0.2273776
#> Pop Std Dev: 0.048747
#> MAD: 0.045125
#> Half 16-84 Quan (1s): 0.045516
#> Half 02-98 Quan (2s): 0.098538
#> Using 4386102 out of 4840000 (90.62%) data points (1304 < xlo & 48632 > xhi)
par(mar=c(3.1,3.1,1.1,1.1))
maghist(stack_clip3$image$imDat, xlim=5)
#> Summary of used sample:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.2027822 -0.0269831 -0.0003092 0.0010894 0.0271422 0.2022182
#> Pop Std Dev: 0.043698
#> MAD: 0.040108
#> Half 16-84 Quan (1s): 0.040507
#> Half 02-98 Quan (2s): 0.08843
#> Using 4378851 out of 4840000 (90.47%) data points (2314 < xlo & 54873 > xhi)
The median stack as a sky RMS near to 0.049 whilst he weighted stack
is deeper at 0.044. The exact ratio is complicated by the clipping
removing effective frames and pixel correlation, but the difference is
still quite notable.
In general median stacking is hard for large images because to do
this the brute force way you need to keep all projected images in memory
for the entire region, so in this case it would be eight 2,200 x 2,200
images, but for large regions (we often stack square degrees) you might
have a few hundred 20,000 x 20,0000. This becomes intractable even on
the biggest memory machines (e.g. that would required TBs of RAM). To
solve this we implement a similar solution to Swarp where projected
frames saved on disk are median stacked in chunks and the contiguous
image is then put back together at the end.
We can make use of the frames we have already dumped to do this with
Rwcs_stack_median:
stack_med2 = propaneStackWarpMed(dirlist = stack_clip$dump_dir,
keyvalues_out = keyvalues_out,
cores = 8,
pattern = glob2rx('*image*fits') #to only stack images
)
#> Time taken: 2.38 seconds
And now we can compare the region with the bad cosmic rays:
plot(stack$image[1000,500,box=500])
plot(stack_clip2$image[1000,500,box=500])
plot(stack_clip3$image[1000,500,box=500])
plot(stack_med2$image[1000,500,box=500])
In general the best quality data will come from using clipped inverse
variance weighted stacks (since it will remove the outliers and keep the
higher \(S/N\)), but if little is known
about the data then a median might be the best practical solution, being
quite robust to artefacts and not requiring data specific clipping
tuning (clip_tol, clip_dilate
etc).
To help automate this patching process we can use the
propanePatch function. Here we use our worst (not clipped)
image and combine with the median stack:
stack_patch = propanePatch(image_inVar = stack$image, image_med = stack_med2$image)
#> Patching 0.14% of threshold pixels and 8.35% of NA pixels
Note we just run with defaults, but there are quite a few options to
determine optimal patching- most importantly the
threshold (where lower mean more patching with the
median combined image). We can also specify just patching
hot or cold pixels (if you are very
certain only certain types of bad pixel exist in the data).
Relatively few very pixels in the stack$image are
replaced with stack_med2$image value using the default arguments
(0.14%). We can see the resulting patched image and check which pixels
were replaced:
plot(stack_patch$image[1000,500,box=500])
plot(stack_patch$patch[1000,500,box=500], qdiff=T)
Now we can see if we see much difference using the already clipped
version:
stack_patch3 = propanePatch(image_inVar = stack_clip3$image, image_med = stack_med2$image)
#> Patching 0.15% of threshold pixels and 8.35% of NA pixels
And we again check the patched plots:
plot(stack_patch3$image[1000,500,box=500])
plot(stack_patch3$patch[1000,500,box=500], qdiff=TRUE)
And finally lets check the differences in the patched images:
plot(stack_patch$image[1000,500,box=500] - stack_patch3$image[1000,500,box=500], qdiff=TRUE)
This would suggest a marginal preference for the clipped and patched
version of the image.