Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):
knitr::opts_chunk$set(eval=FALSE)
Load some packages.
library(ProFound) library(Rfits) library(Rwcs)
First we load some images.
im_r = Rfits_read_image(system.file("extdata", 'MultiBand/r.fits', package="ProFound"),ext=2)$imDat im_i = Rfits_read_image(system.file("extdata", 'MultiBand/i.fits', package="ProFound"),ext=2)$imDat im_Z = Rfits_read_image(system.file("extdata", 'MultiBand/Z.fits', package="ProFound"),ext=2)$imDat im_Y = Rfits_read_image(system.file("extdata", 'MultiBand/Y.fits', package="ProFound"),ext=2)$imDat
Check them out.
magimage(im_r) magimage(im_i) magimage(im_Z) magimage(im_Y)
We run ProFound in different modes in the optical and NIR. In the optical we want to detect right into the outskirts of faint extended sources, and in the NIR we want to avoid finding spurious objects in the very noisy regions we occassionally see (this is for WAVES image analysis, this particular frame is okay).
pro_r = profoundProFound(im_r, skycut=2, roughpedestal=T, grid=50, boxiters=4, verbose=F, SBdilate=2, box=200, stats=FALSE) pro_i = profoundProFound(im_i, skycut=2, roughpedestal=T, grid=50, boxiters=4, verbose=F, SBdilate=2, box=200, stats=FALSE) pro_Z = profoundProFound(im_Z, skycut=5, sigma=5, roughpedestal=T, grid=50, boxiters=4, verbose=F, SBdilate=2, box=200, magzero=30, stats=FALSE) pro_Y = profoundProFound(im_Y, skycut=5, sigma=5, roughpedestal=T, grid=50, boxiters=4, verbose=F, SBdilate=2, box=200, magzero=30, stats=FALSE)
We create a new sky map using a finer background to find the sharp undesired features.
sky_r = profoundMakeSkyGrid(im_r, objects=pro_r$objects_redo, doclip=T, grid=50, box=50, boxiters=4) sky_i = profoundMakeSkyGrid(im_i, objects=pro_i$objects_redo, doclip=T, grid=50, box=50, boxiters=4) sky_Z = profoundMakeSkyGrid(im_Z, objects=pro_Z$objects_redo, doclip=T, grid=50, box=50, boxiters=4) sky_Y = profoundMakeSkyGrid(im_Y, objects=pro_Y$objects_redo, doclip=T, grid=50, box=50, boxiters=4)
We mask pixels out based on global statistics. For the optical we want to avoid the detector edges (so mask based on sky), whereas in the NIR we want to avoid the really noisy regions (so mask based on sky RMS).
mask_r = sky_r$sky > (median(sky_r$sky,na.rm=TRUE) + 0.5*(median(sky_r$sky,,na.rm=TRUE) - quantile(sky_r$sky,0.01,na.rm=TRUE))) mask_i = sky_i$sky > (median(sky_i$sky,na.rm=TRUE) + 1*(median(sky_i$sky,na.rm=TRUE) - quantile(sky_i$sky,0.01,na.rm=TRUE))) mask_Z = sky_Z$skyRMS > (median(sky_Z$skyRMS,na.rm=TRUE) + 1*(median(sky_Z$skyRMS,na.rm=TRUE) - quantile(sky_Z$skyRMS,0.01,na.rm=TRUE))) mask_Y = sky_Y$skyRMS > (median(sky_Y$skyRMS,na.rm=TRUE) + 1*(median(sky_Y$skyRMS,na.rm=TRUE) - quantile(sky_Y$skyRMS,0.01,na.rm=TRUE)))
We unmask based on a highly blurred r band image.
unmask = profoundImBlur(im_r,20) > 1e-11 mask_r[unmask] = 0 mask_i[unmask] = 0 mask_Z[unmask] = 0 mask_Y[unmask] = 0
magimage(pro_r$image); magimage(mask_r, magmap=F, col=c(NA,hsv(alpha=0.3)),add=TRUE) magimage(pro_i$image); magimage(mask_i, magmap=F, col=c(NA,hsv(alpha=0.3)),add=TRUE) magimage(pro_Z$image); magimage(mask_Z, magmap=F, col=c(NA,hsv(alpha=0.3)),add=TRUE) magimage(pro_Y$image); magimage(mask_Y, magmap=F, col=c(NA,hsv(alpha=0.3)),add=TRUE)
Now we can do some stacking. We try rZ, riZ and riZY variants.
stack_rZ = profoundMakeStack(image_list = list(pro_r$image, pro_Z$image), sky_list = list(pro_r$sky, pro_Z$sky), skyRMS_list = list(pro_r$skyRMS, pro_Z$skyRMS), mask_list = list(mask_r, mask_Z), magzero_in = c(0,30) )
stack_riZ = profoundMakeStack(image_list = list(pro_r$image, pro_i$image, pro_Z$image), sky_list = list(pro_r$sky, pro_i$sky, pro_Z$sky), skyRMS_list = list(pro_r$skyRMS, pro_i$skyRMS, pro_Z$skyRMS), mask_list = list(mask_r, mask_i, mask_Z), magzero_in = c(0,0,30) )
stack_riZY = profoundMakeStack(image_list = list(pro_r$image, pro_i$image, pro_Z$image, pro_Y$image), sky_list = list(pro_r$sky, pro_i$sky, pro_Z$sky, pro_Y$sky), skyRMS_list = list(pro_r$skyRMS, pro_i$skyRMS, pro_Z$skyRMS, pro_Y$skyRMS), mask_list = list(mask_r, mask_i, mask_Z, mask_Y), magzero_in = c(0,0,30,30) )
Check out the different stacks. Overall, riZY looks best.
magimage(stack_rZ$image) magimage(stack_riZ$image) magimage(stack_riZY$image)
Run ProFound on the different detection images.
pro_detect_rZ = profoundProFound(stack_rZ$image, skycut=2, sky=0, grid=50, boxiters=4, SBdilate=2, tolerance=2, sigma=0.8) pro_detect_riZ = profoundProFound(stack_riZ$image, skycut=2, sky=0, grid=50, boxiters=4, SBdilate=2, tolerance=2, sigma=0.8) pro_detect_riZY = profoundProFound(stack_riZY$image, skycut=2, sky=0, grid=50, boxiters=4, SBdilate=2, tolerance=2, sigma=0.8)
Check the summary plots.
plot(pro_detect_rZ) plot(pro_detect_riZ) plot(pro_detect_riZY)
Zoom in to an interesting region showing a variety of difficult detection issues. Again, riZY look to do the best overall.
profoundSegimPlot(pro_detect_rZ$image[1000+ -500:499,700+ -500:499], segim=pro_detect_rZ$segim[1000+ -500:499,700+ -500:499]) profoundSegimPlot(pro_detect_riZ$image[1000+ -500:499,700+ -500:499], segim=pro_detect_riZ$segim[1000+ -500:499,700+ -500:499]) profoundSegimPlot(pro_detect_riZY$image[1000+ -500:499,700+ -500:499], segim=pro_detect_riZY$segim[1000+ -500:499,700+ -500:499])
pro_r_detect = profoundProFound(im_r, segim=pro_detect_rZ$segim, sky=pro_r$sky, skyRMS=pro_r$skyRMS, static_photom=TRUE) pro_Z_detect = profoundProFound(im_Z, segim=pro_detect_rZ$segim, sky=pro_Z$sky, skyRMS=pro_Z$skyRMS, static_photom=TRUE, magzero=30)
segstats = cbind(pro_r_detect$segstats, col=pro_r_detect$segstats$mag-pro_Z_detect$segstats$mag)
groups = profoundAutoMerge(pro_detect_rZ$segim_orig, segstats = segstats, spur_lim=4e-3, col_lim = c(0,0.8)) segim_fix = profoundSegimKeep(pro_detect_rZ$segim, segID_merge=groups$segID)
profoundSegimPlot(pro_detect_riZY$image, segim=segim_fix)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.