ProPane Stacking: Let’s do the Space Warp Again!

Aaron Robotham

2026-04-16

First we load the packages we will need for this vignette:

library(ProPane)
library(Rfits)
library(Rwcs)
library(foreach)
library(magicaxis)
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

Set the number of cores to use:

cores = 8L
multitype = 'cluster'

The ProPane package comes with key utilities for warping between different WCS systems: propaneWarp (for warping individual frames once). ProPane also contains the various functions for creating large stacks of many warped frames (which is of class ProPane, which is roughly meant to suggest the idea of many panes of glass being stacked together).

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 = cores, multitype = multitype,
                   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: 10.57 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 
#> 0.08346322314049587 0.01885909090909091 0.00711466942148760 0.00850681818181818 
#>                   4                   5                   6                   7 
#> 0.01295826446280992 0.00806859504132231 0.01360061983471074 0.01583677685950413 
#>                   8 
#> 0.83159194214876031

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 
#> -0.196649352849538 -0.026858259933775 -0.000384779622841  0.001026429638022 
#>            3rd Qu.               Max. 
#>  0.026834725731099  0.199820488220712
#> Pop Std Dev: 0.042005
#> MAD: 0.039786
#> Half 16-84 Quan (1s): 0.040041
#> Half 02-98 Quan (2s): 0.084312
#> Using 3974591 out of 4024905 (98.75%) data points (9 < xlo & 50305 > 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.0458381112992776

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 cast 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 = cores, multitype = multitype,
                   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: 7.366 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 = cores, multitype = multitype,
                   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: 8.056 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 
#> -1.36430564719320 -0.18659601673452 -0.01085886655576  0.00436604341461 
#>           3rd Qu.              Max. 
#>  0.17300958460594  1.34320361445168
#> Pop Std Dev: 0.30086
#> MAD: 0.26642
#> Half 16-84 Quan (1s): 0.27091
#> Half 02-98 Quan (2s): 0.62411
#> Using 393248 out of 640000 (61.44%) data points (258 < xlo & 7849 > 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 
#> -2.30884842203178 -0.30833700359419 -0.00216663944002  0.01373513601945 
#>           3rd Qu.              Max. 
#>  0.31125336617209  2.30966847980062
#> Pop Std Dev: 0.499
#> MAD: 0.45901
#> Half 16-84 Quan (1s): 0.46238
#> Half 02-98 Quan (2s): 1.0099
#> 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.410452638014277
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.289893959921438

The default is to project from the higher resolution grid, so in this case without specifying 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 = cores, multitype = multitype,
                   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//RtmpJHi3xC
#> 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
#> Called from: propaneStackWarpInVar(image_list = image_list, inVar_list = 0.000468, 
#>     exp_list = 10, cores = cores, multitype = multitype, keyvalues_out = keyvalues_out, 
#>     magzero_in = 30, magzero_out = 23.9, keep_extreme_pix = TRUE, 
#>     keepcrop = FALSE, doclip = TRUE, dump_frames = TRUE, return_all = TRUE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#595: post_stack_DT = data.table(cold = as.integer(post_stack_cold_id), 
#>     hot = as.integer(post_stack_hot_id), ID = 1:length(post_stack_cold_id), 
#>     key = c("cold", "hot"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#600: rm(post_stack_cold_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#601: rm(post_stack_hot_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#604: if (Nbatch == Nim) {
#>     message("Restacking without clipped cold/hot pixels")
#>     if (is.null(pre_stack_inVar_list)) {
#>         message("Stacking Images ", seq_start, " to ", seq_end, 
#>             " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#>     else {
#>         message("Stacking Images and InVar ", seq_start, " to ", 
#>             seq_end, " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                 post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> } else {
#>     message("Reprojecting and restacking without clipped cold/hot pixels")
#>     for (seq_start in seq_process) {
#>         seq_end = min(seq_start + Nbatch - 1L, Nim)
#>         Nbatch_sub = length(seq_start:seq_end)
#>         pre_stack_image_list = NULL
#>         pre_stack_inVar_list = NULL
#>         pre_stack_exp_list = NULL
#>         pre_stack_weight_list = NULL
#>         if (dump_frames) {
#>             message("Loading Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end) %dopar% 
#>                 {
#>                   Rfits_read_image(paste0(dump_dir, "/image_warp_", 
#>                     i, ".fits"))
#>                 }
#>         }
#>         else {
#>             message("Projecting Images ", seq_start, " to ", 
#>                 seq_end, " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end, 
#>                 .noexport = c("post_stack_image", "post_stack_weight", 
#>                   "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                 {
#>                   needs_read = !is.null(mask_list) || !is.null(inVar_list)
#>                   if (inherits(image_list[[i]], "Rfits_pointer") && 
#>                     needs_read) {
#>                     temp_image = image_list[[i]][, ]
#>                   }
#>                   else {
#>                     temp_image = image_list[[i]]
#>                   }
#>                   if (!is.null(mask_list)) {
#>                     if (inherits(mask_list[[i]], "Rfits_pointer")) {
#>                       mask_list[[i]]$header = FALSE
#>                       temp_image$imDat[mask_list[[i]][, ] != 
#>                         0] = NA
#>                     }
#>                     else if (inherits(mask_list[[i]], "Rfits_image")) {
#>                       temp_image$imDat[mask_list[[i]]$imDat != 
#>                         0] = NA
#>                     }
#>                     else {
#>                       temp_image$imDat[mask_list != 0] = NA
#>                     }
#>                   }
#>                   if (!is.null(inVar_list)) {
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar = inVar_list[[i]][, ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar = inVar_list[[i]]
#>                     }
#>                     if (is.matrix(temp_inVar)) {
#>                       temp_image$imDat[temp_inVar == 0] = NA
#>                     }
#>                     rm(temp_inVar)
#>                   }
#>                   suppressWarnings({
#>                     temp_warp = propaneWarp(image_in = temp_image, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = TRUE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield_return = TRUE, cores = cores_warp, 
#>                       ...)
#>                     if (zero_point_scale[i] != 1) {
#>                       temp_warp$imDat = temp_warp$imDat * zero_point_scale[i]
#>                     }
#>                   })
#>                   return(temp_warp)
#>                 }
#>         }
#>         if (!is.null(inVar_list)) {
#>             if (dump_frames) {
#>                 message("Loading Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end) %dopar% 
#>                   {
#>                     Rfits_read_image(paste0(dump_dir, "/inVar_warp_", 
#>                       i, ".fits"))
#>                   }
#>             }
#>             else {
#>                 message("Projecting Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end, 
#>                   .noexport = c("post_stack_image", "post_stack_weight", 
#>                     "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                   {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_inVar = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_inVar = image_list[[i]]
#>                     }
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar$imDat[] = inVar_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar$imDat[] = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar$imDat[] = inVar_list[[i]]
#>                     }
#>                     suppressMessages({
#>                       temp_warp = propaneWarp(image_in = temp_inVar, 
#>                         keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                         doscale = FALSE, dotightcrop = TRUE, 
#>                         keepcrop = keepcrop, warpfield = pre_stack_image_list[[i - 
#>                           seq_start + 1L]]$warpfield, cores = cores_warp, 
#>                         ...)
#>                       if (zero_point_scale[i] != 1) {
#>                         temp_warp$imDat = temp_warp$imDat/(zero_point_scale[i]^2)
#>                       }
#>                       temp_warp$imDat = temp_warp$imDat * (Rwcs_pixscale(temp_inVar$keyvalues)^4/Rwcs_pixscale(keyvalues_out)^4)
#>                       return(temp_warp)
#>                     })
#>                   }
#>             }
#>         }
#>         if (any(weight_image)) {
#>             if (dump_frames) {
#>                 message("Loading Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             else {
#>                 message("Projecting Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             pre_stack_weight_list = foreach(i = seq_start:seq_end, 
#>                 .packages = c("Rwcs", "Rfits"), export = c("image_list", 
#>                   "dump_dir", "weight_image")) %dopar% {
#>                 if (weight_image[i]) {
#>                   if (dump_frames) {
#>                     Rfits_read_image(paste0(dump_dir, "/weight_warp_", 
#>                       i, ".fits"))
#>                   }
#>                   else {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_weight = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_weight = image_list[[i]]
#>                     }
#>                     if (inherits(weight_list[[i]], "Rfits_pointer")) {
#>                       weight_list[[i]]$header = FALSE
#>                       temp_weight$imDat[] = weight_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(weight_list[[i]], "Rfits_image")) {
#>                       temp_weight$imDat[] = weight_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_weight$imDat[] = weight_list[[i]]
#>                     }
#>                     return(propaneWarp(image_in = temp_weight, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = FALSE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield = pre_stack_image_list[[i - seq_start + 
#>                         1L]]$warpfield, cores = cores_warp, ...))
#>                   }
#>                 }
#>                 else {
#>                   return(weight_list[[i]])
#>                 }
#>             }
#>         }
#>         else {
#>             pre_stack_weight_list = weight_list
#>         }
#>         if (is.null(pre_stack_inVar_list)) {
#>             message("Stacking Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_cpp(post_image = post_stack_image, 
#>                   post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         else {
#>             message("Stacking Images and InVar ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                   post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                   pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         if (keep_extreme_pix) {
#>             message("Calculating Extreme Pixels ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             temp_mask_clip = matrix(0L, dim_im[1], dim_im[2])
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mode(temp_mask_clip) = "logical"
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#605: message("Restacking without clipped cold/hot pixels")
#> Restacking without clipped cold/hot pixels
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#606: if (is.null(pre_stack_inVar_list)) {
#>     message("Stacking Images ", seq_start, " to ", seq_end, " of ", 
#>         Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>             offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> } else {
#>     message("Stacking Images and InVar ", seq_start, " to ", 
#>         seq_end, " of ", Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_inVar_cpp(post_image = post_stack_image, 
#>             post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>             pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#658: message("Stacking Images and InVar ", seq_start, " to ", seq_end, 
#>     " of ", Nim)
#> Stacking Images and InVar 1 to 8 of 8
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#659: for (i in 1:Nim) {
#>     if (i > 1) {
#>         temp_mask_clip[] = 0L
#>     }
#>     temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>         i, ID]))] = 1L
#>     if (clip_dilate > 0) {
#>         temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>     }
#>     mask_clip = mask_clip + temp_mask_clip
#>     mode(temp_mask_clip) = "logical"
#>     if (dump_frames) {
#>         Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>             i, ".fits"))
#>     }
#>     if (weight_image[i]) {
#>         pre_weight = pre_stack_weight_list[[i]]$imDat
#>     }
#>     else {
#>         pre_weight = weight_list[[i]]
#>     }
#>     .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>         post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>         pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>         offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>             "YCUTLO")]), post_mask = temp_mask_clip)
#>     if (keep_extreme_pix) {
#>         xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>         ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>         new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>         post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#672: Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>     i, ".fits"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#993: if (return_all == FALSE) {
#>     pre_stack_exp_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#997: weight_sel = (post_stack_weight != 0L)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#999: if (is.null(pre_stack_inVar_list)) {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_weight[weight_sel]
#> } else {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#>     post_stack_inVar[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1002: post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1003: post_stack_inVar[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1006: if (keep_extreme_pix) {
#>     post_stack_cold[!weight_sel] = NA
#>     post_stack_hot[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1007: post_stack_cold[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1008: post_stack_hot[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1011: post_stack_image[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1014: if (return_all == FALSE) {
#>     pre_stack_image_list = NULL
#>     pre_stack_inVar_list = NULL
#>     pre_stack_exp_list = NULL
#>     pre_stack_weight_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1021: keyvalues_out$EXTNAME = "image"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1022: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1023: keyvalues_out$R_VER = R.version$version.string
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1024: keyvalues_out$PANE_VER = as.character(packageVersion("ProPane"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1025: keyvalues_out$RWCS_VER = as.character(packageVersion("Rwcs"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1027: image_out = Rfits_create_image(post_stack_image, keyvalues = keyvalues_out, 
#>     keypass = FALSE, history = "Stacked with propaneStackWarpInVar")
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1032: keyvalues_out$EXTNAME = "weight"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1033: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1034: weight_out = Rfits_create_image(post_stack_weight, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1038: if (!is.null(post_stack_inVar)) {
#>     keyvalues_out$EXTNAME = "inVar"
#>     keyvalues_out$MAGZERO = magzero_out
#>     inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     inVar_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1039: keyvalues_out$EXTNAME = "inVar"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1040: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1041: inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1048: if (!is.null(post_stack_exp)) {
#>     keyvalues_out$EXTNAME = "exp"
#>     keyvalues_out$MAGZERO = NULL
#>     exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     exp_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1049: keyvalues_out$EXTNAME = "exp"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1050: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1051: exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1058: if (keep_extreme_pix) {
#>     post_stack_cold[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "cold"
#>     keyvalues_out$MAGZERO = magzero_out
#>     cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#>     post_stack_hot[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "hot"
#>     hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     cold_out = NULL
#>     hot_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1059: post_stack_cold[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1061: keyvalues_out$EXTNAME = "cold"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1062: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1063: cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1067: post_stack_hot[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1069: keyvalues_out$EXTNAME = "hot"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1070: hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1078: if (!is.null(mask_clip)) {
#>     mask_clip[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "clip"
#>     keyvalues_out$MAGZERO = NULL
#>     mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1079: mask_clip[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1081: keyvalues_out$EXTNAME = "clip"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1082: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1083: mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1088: time_taken = proc.time()[3] - timestart
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1089: message("Time taken: ", signif(time_taken, 4), " seconds")
#> Time taken: 17.62 seconds
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1091: if (return_all) {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>         exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> } else {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1092: output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>     exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>     image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>     exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>     which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>     dump_dir = dump_dir)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1120: class(output) = "ProPane"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1121: stopImplicitCluster()
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1122: return(invisible(output))

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 = cores, multitype = multitype,
                   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
#> Called from: propaneStackWarpInVar(image_list = image_list, inVar_list = 0.000468, 
#>     exp_list = 10, cores = cores, multitype = multitype, keyvalues_out = keyvalues_out, 
#>     magzero_in = 30, magzero_out = 23.9, keep_extreme_pix = TRUE, 
#>     keepcrop = FALSE, doclip = TRUE, clip_tol = 10, dump_frames = FALSE, 
#>     return_all = TRUE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#595: post_stack_DT = data.table(cold = as.integer(post_stack_cold_id), 
#>     hot = as.integer(post_stack_hot_id), ID = 1:length(post_stack_cold_id), 
#>     key = c("cold", "hot"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#600: rm(post_stack_cold_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#601: rm(post_stack_hot_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#604: if (Nbatch == Nim) {
#>     message("Restacking without clipped cold/hot pixels")
#>     if (is.null(pre_stack_inVar_list)) {
#>         message("Stacking Images ", seq_start, " to ", seq_end, 
#>             " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#>     else {
#>         message("Stacking Images and InVar ", seq_start, " to ", 
#>             seq_end, " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                 post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> } else {
#>     message("Reprojecting and restacking without clipped cold/hot pixels")
#>     for (seq_start in seq_process) {
#>         seq_end = min(seq_start + Nbatch - 1L, Nim)
#>         Nbatch_sub = length(seq_start:seq_end)
#>         pre_stack_image_list = NULL
#>         pre_stack_inVar_list = NULL
#>         pre_stack_exp_list = NULL
#>         pre_stack_weight_list = NULL
#>         if (dump_frames) {
#>             message("Loading Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end) %dopar% 
#>                 {
#>                   Rfits_read_image(paste0(dump_dir, "/image_warp_", 
#>                     i, ".fits"))
#>                 }
#>         }
#>         else {
#>             message("Projecting Images ", seq_start, " to ", 
#>                 seq_end, " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end, 
#>                 .noexport = c("post_stack_image", "post_stack_weight", 
#>                   "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                 {
#>                   needs_read = !is.null(mask_list) || !is.null(inVar_list)
#>                   if (inherits(image_list[[i]], "Rfits_pointer") && 
#>                     needs_read) {
#>                     temp_image = image_list[[i]][, ]
#>                   }
#>                   else {
#>                     temp_image = image_list[[i]]
#>                   }
#>                   if (!is.null(mask_list)) {
#>                     if (inherits(mask_list[[i]], "Rfits_pointer")) {
#>                       mask_list[[i]]$header = FALSE
#>                       temp_image$imDat[mask_list[[i]][, ] != 
#>                         0] = NA
#>                     }
#>                     else if (inherits(mask_list[[i]], "Rfits_image")) {
#>                       temp_image$imDat[mask_list[[i]]$imDat != 
#>                         0] = NA
#>                     }
#>                     else {
#>                       temp_image$imDat[mask_list != 0] = NA
#>                     }
#>                   }
#>                   if (!is.null(inVar_list)) {
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar = inVar_list[[i]][, ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar = inVar_list[[i]]
#>                     }
#>                     if (is.matrix(temp_inVar)) {
#>                       temp_image$imDat[temp_inVar == 0] = NA
#>                     }
#>                     rm(temp_inVar)
#>                   }
#>                   suppressWarnings({
#>                     temp_warp = propaneWarp(image_in = temp_image, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = TRUE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield_return = TRUE, cores = cores_warp, 
#>                       ...)
#>                     if (zero_point_scale[i] != 1) {
#>                       temp_warp$imDat = temp_warp$imDat * zero_point_scale[i]
#>                     }
#>                   })
#>                   return(temp_warp)
#>                 }
#>         }
#>         if (!is.null(inVar_list)) {
#>             if (dump_frames) {
#>                 message("Loading Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end) %dopar% 
#>                   {
#>                     Rfits_read_image(paste0(dump_dir, "/inVar_warp_", 
#>                       i, ".fits"))
#>                   }
#>             }
#>             else {
#>                 message("Projecting Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end, 
#>                   .noexport = c("post_stack_image", "post_stack_weight", 
#>                     "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                   {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_inVar = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_inVar = image_list[[i]]
#>                     }
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar$imDat[] = inVar_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar$imDat[] = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar$imDat[] = inVar_list[[i]]
#>                     }
#>                     suppressMessages({
#>                       temp_warp = propaneWarp(image_in = temp_inVar, 
#>                         keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                         doscale = FALSE, dotightcrop = TRUE, 
#>                         keepcrop = keepcrop, warpfield = pre_stack_image_list[[i - 
#>                           seq_start + 1L]]$warpfield, cores = cores_warp, 
#>                         ...)
#>                       if (zero_point_scale[i] != 1) {
#>                         temp_warp$imDat = temp_warp$imDat/(zero_point_scale[i]^2)
#>                       }
#>                       temp_warp$imDat = temp_warp$imDat * (Rwcs_pixscale(temp_inVar$keyvalues)^4/Rwcs_pixscale(keyvalues_out)^4)
#>                       return(temp_warp)
#>                     })
#>                   }
#>             }
#>         }
#>         if (any(weight_image)) {
#>             if (dump_frames) {
#>                 message("Loading Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             else {
#>                 message("Projecting Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             pre_stack_weight_list = foreach(i = seq_start:seq_end, 
#>                 .packages = c("Rwcs", "Rfits"), export = c("image_list", 
#>                   "dump_dir", "weight_image")) %dopar% {
#>                 if (weight_image[i]) {
#>                   if (dump_frames) {
#>                     Rfits_read_image(paste0(dump_dir, "/weight_warp_", 
#>                       i, ".fits"))
#>                   }
#>                   else {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_weight = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_weight = image_list[[i]]
#>                     }
#>                     if (inherits(weight_list[[i]], "Rfits_pointer")) {
#>                       weight_list[[i]]$header = FALSE
#>                       temp_weight$imDat[] = weight_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(weight_list[[i]], "Rfits_image")) {
#>                       temp_weight$imDat[] = weight_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_weight$imDat[] = weight_list[[i]]
#>                     }
#>                     return(propaneWarp(image_in = temp_weight, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = FALSE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield = pre_stack_image_list[[i - seq_start + 
#>                         1L]]$warpfield, cores = cores_warp, ...))
#>                   }
#>                 }
#>                 else {
#>                   return(weight_list[[i]])
#>                 }
#>             }
#>         }
#>         else {
#>             pre_stack_weight_list = weight_list
#>         }
#>         if (is.null(pre_stack_inVar_list)) {
#>             message("Stacking Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_cpp(post_image = post_stack_image, 
#>                   post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         else {
#>             message("Stacking Images and InVar ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                   post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                   pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         if (keep_extreme_pix) {
#>             message("Calculating Extreme Pixels ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             temp_mask_clip = matrix(0L, dim_im[1], dim_im[2])
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mode(temp_mask_clip) = "logical"
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#605: message("Restacking without clipped cold/hot pixels")
#> Restacking without clipped cold/hot pixels
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#606: if (is.null(pre_stack_inVar_list)) {
#>     message("Stacking Images ", seq_start, " to ", seq_end, " of ", 
#>         Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>             offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> } else {
#>     message("Stacking Images and InVar ", seq_start, " to ", 
#>         seq_end, " of ", Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_inVar_cpp(post_image = post_stack_image, 
#>             post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>             pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#658: message("Stacking Images and InVar ", seq_start, " to ", seq_end, 
#>     " of ", Nim)
#> Stacking Images and InVar 1 to 8 of 8
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#659: for (i in 1:Nim) {
#>     if (i > 1) {
#>         temp_mask_clip[] = 0L
#>     }
#>     temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>         i, ID]))] = 1L
#>     if (clip_dilate > 0) {
#>         temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>     }
#>     mask_clip = mask_clip + temp_mask_clip
#>     mode(temp_mask_clip) = "logical"
#>     if (dump_frames) {
#>         Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>             i, ".fits"))
#>     }
#>     if (weight_image[i]) {
#>         pre_weight = pre_stack_weight_list[[i]]$imDat
#>     }
#>     else {
#>         pre_weight = weight_list[[i]]
#>     }
#>     .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>         post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>         pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>         offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>             "YCUTLO")]), post_mask = temp_mask_clip)
#>     if (keep_extreme_pix) {
#>         xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>         ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>         new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>         post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#993: if (return_all == FALSE) {
#>     pre_stack_exp_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#997: weight_sel = (post_stack_weight != 0L)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#999: if (is.null(pre_stack_inVar_list)) {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_weight[weight_sel]
#> } else {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#>     post_stack_inVar[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1002: post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1003: post_stack_inVar[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1006: if (keep_extreme_pix) {
#>     post_stack_cold[!weight_sel] = NA
#>     post_stack_hot[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1007: post_stack_cold[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1008: post_stack_hot[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1011: post_stack_image[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1014: if (return_all == FALSE) {
#>     pre_stack_image_list = NULL
#>     pre_stack_inVar_list = NULL
#>     pre_stack_exp_list = NULL
#>     pre_stack_weight_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1021: keyvalues_out$EXTNAME = "image"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1022: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1023: keyvalues_out$R_VER = R.version$version.string
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1024: keyvalues_out$PANE_VER = as.character(packageVersion("ProPane"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1025: keyvalues_out$RWCS_VER = as.character(packageVersion("Rwcs"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1027: image_out = Rfits_create_image(post_stack_image, keyvalues = keyvalues_out, 
#>     keypass = FALSE, history = "Stacked with propaneStackWarpInVar")
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1032: keyvalues_out$EXTNAME = "weight"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1033: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1034: weight_out = Rfits_create_image(post_stack_weight, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1038: if (!is.null(post_stack_inVar)) {
#>     keyvalues_out$EXTNAME = "inVar"
#>     keyvalues_out$MAGZERO = magzero_out
#>     inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     inVar_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1039: keyvalues_out$EXTNAME = "inVar"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1040: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1041: inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1048: if (!is.null(post_stack_exp)) {
#>     keyvalues_out$EXTNAME = "exp"
#>     keyvalues_out$MAGZERO = NULL
#>     exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     exp_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1049: keyvalues_out$EXTNAME = "exp"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1050: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1051: exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1058: if (keep_extreme_pix) {
#>     post_stack_cold[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "cold"
#>     keyvalues_out$MAGZERO = magzero_out
#>     cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#>     post_stack_hot[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "hot"
#>     hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     cold_out = NULL
#>     hot_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1059: post_stack_cold[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1061: keyvalues_out$EXTNAME = "cold"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1062: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1063: cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1067: post_stack_hot[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1069: keyvalues_out$EXTNAME = "hot"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1070: hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1078: if (!is.null(mask_clip)) {
#>     mask_clip[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "clip"
#>     keyvalues_out$MAGZERO = NULL
#>     mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1079: mask_clip[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1081: keyvalues_out$EXTNAME = "clip"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1082: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1083: mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1088: time_taken = proc.time()[3] - timestart
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1089: message("Time taken: ", signif(time_taken, 4), " seconds")
#> Time taken: 15.85 seconds
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1091: if (return_all) {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>         exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> } else {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1092: output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>     exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>     image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>     exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>     which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>     dump_dir = dump_dir)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1120: class(output) = "ProPane"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1121: stopImplicitCluster()
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1122: return(invisible(output))
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 = cores, multitype = multitype,
                   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
#> Called from: propaneStackWarpInVar(image_list = image_list, inVar_list = 0.000468, 
#>     exp_list = 10, cores = cores, multitype = multitype, keyvalues_out = keyvalues_out, 
#>     magzero_in = 30, magzero_out = 23.9, keep_extreme_pix = TRUE, 
#>     keepcrop = FALSE, doclip = TRUE, clip_tol = 50, clip_dilate = 5, 
#>     return_all = TRUE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#595: post_stack_DT = data.table(cold = as.integer(post_stack_cold_id), 
#>     hot = as.integer(post_stack_hot_id), ID = 1:length(post_stack_cold_id), 
#>     key = c("cold", "hot"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#600: rm(post_stack_cold_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#601: rm(post_stack_hot_id)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#604: if (Nbatch == Nim) {
#>     message("Restacking without clipped cold/hot pixels")
#>     if (is.null(pre_stack_inVar_list)) {
#>         message("Stacking Images ", seq_start, " to ", seq_end, 
#>             " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#>     else {
#>         message("Stacking Images and InVar ", seq_start, " to ", 
#>             seq_end, " of ", Nim)
#>         for (i in 1:Nim) {
#>             if (i > 1) {
#>                 temp_mask_clip[] = 0L
#>             }
#>             temp_mask_clip[unique(c(post_stack_DT[cold == i, 
#>                 ID], post_stack_DT[hot == i, ID]))] = 1L
#>             if (clip_dilate > 0) {
#>                 temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>             }
#>             mask_clip = mask_clip + temp_mask_clip
#>             mode(temp_mask_clip) = "logical"
#>             if (dump_frames) {
#>                 Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                   "/mask_warp_", i, ".fits"))
#>             }
#>             if (weight_image[i]) {
#>                 pre_weight = pre_stack_weight_list[[i]]$imDat
#>             }
#>             else {
#>                 pre_weight = weight_list[[i]]
#>             }
#>             .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                 post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                 pre_image = pre_stack_image_list[[i]]$imDat, 
#>                 pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                 pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                   "YCUTLO")]), post_mask = temp_mask_clip)
#>             if (keep_extreme_pix) {
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> } else {
#>     message("Reprojecting and restacking without clipped cold/hot pixels")
#>     for (seq_start in seq_process) {
#>         seq_end = min(seq_start + Nbatch - 1L, Nim)
#>         Nbatch_sub = length(seq_start:seq_end)
#>         pre_stack_image_list = NULL
#>         pre_stack_inVar_list = NULL
#>         pre_stack_exp_list = NULL
#>         pre_stack_weight_list = NULL
#>         if (dump_frames) {
#>             message("Loading Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end) %dopar% 
#>                 {
#>                   Rfits_read_image(paste0(dump_dir, "/image_warp_", 
#>                     i, ".fits"))
#>                 }
#>         }
#>         else {
#>             message("Projecting Images ", seq_start, " to ", 
#>                 seq_end, " of ", Nim)
#>             pre_stack_image_list = foreach(i = seq_start:seq_end, 
#>                 .noexport = c("post_stack_image", "post_stack_weight", 
#>                   "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                 {
#>                   needs_read = !is.null(mask_list) || !is.null(inVar_list)
#>                   if (inherits(image_list[[i]], "Rfits_pointer") && 
#>                     needs_read) {
#>                     temp_image = image_list[[i]][, ]
#>                   }
#>                   else {
#>                     temp_image = image_list[[i]]
#>                   }
#>                   if (!is.null(mask_list)) {
#>                     if (inherits(mask_list[[i]], "Rfits_pointer")) {
#>                       mask_list[[i]]$header = FALSE
#>                       temp_image$imDat[mask_list[[i]][, ] != 
#>                         0] = NA
#>                     }
#>                     else if (inherits(mask_list[[i]], "Rfits_image")) {
#>                       temp_image$imDat[mask_list[[i]]$imDat != 
#>                         0] = NA
#>                     }
#>                     else {
#>                       temp_image$imDat[mask_list != 0] = NA
#>                     }
#>                   }
#>                   if (!is.null(inVar_list)) {
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar = inVar_list[[i]][, ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar = inVar_list[[i]]
#>                     }
#>                     if (is.matrix(temp_inVar)) {
#>                       temp_image$imDat[temp_inVar == 0] = NA
#>                     }
#>                     rm(temp_inVar)
#>                   }
#>                   suppressWarnings({
#>                     temp_warp = propaneWarp(image_in = temp_image, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = TRUE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield_return = TRUE, cores = cores_warp, 
#>                       ...)
#>                     if (zero_point_scale[i] != 1) {
#>                       temp_warp$imDat = temp_warp$imDat * zero_point_scale[i]
#>                     }
#>                   })
#>                   return(temp_warp)
#>                 }
#>         }
#>         if (!is.null(inVar_list)) {
#>             if (dump_frames) {
#>                 message("Loading Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end) %dopar% 
#>                   {
#>                     Rfits_read_image(paste0(dump_dir, "/inVar_warp_", 
#>                       i, ".fits"))
#>                   }
#>             }
#>             else {
#>                 message("Projecting Inverse Variance ", seq_start, 
#>                   " to ", seq_end, " of ", Nim)
#>                 pre_stack_inVar_list = foreach(i = seq_start:seq_end, 
#>                   .noexport = c("post_stack_image", "post_stack_weight", 
#>                     "post_stack_inVar", "post_stack_exp")) %dopar% 
#>                   {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_inVar = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_inVar = image_list[[i]]
#>                     }
#>                     if (inherits(inVar_list[[i]], "Rfits_pointer")) {
#>                       inVar_list[[i]]$header = FALSE
#>                       temp_inVar$imDat[] = inVar_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(inVar_list[[i]], "Rfits_image")) {
#>                       temp_inVar$imDat[] = inVar_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_inVar$imDat[] = inVar_list[[i]]
#>                     }
#>                     suppressMessages({
#>                       temp_warp = propaneWarp(image_in = temp_inVar, 
#>                         keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                         doscale = FALSE, dotightcrop = TRUE, 
#>                         keepcrop = keepcrop, warpfield = pre_stack_image_list[[i - 
#>                           seq_start + 1L]]$warpfield, cores = cores_warp, 
#>                         ...)
#>                       if (zero_point_scale[i] != 1) {
#>                         temp_warp$imDat = temp_warp$imDat/(zero_point_scale[i]^2)
#>                       }
#>                       temp_warp$imDat = temp_warp$imDat * (Rwcs_pixscale(temp_inVar$keyvalues)^4/Rwcs_pixscale(keyvalues_out)^4)
#>                       return(temp_warp)
#>                     })
#>                   }
#>             }
#>         }
#>         if (any(weight_image)) {
#>             if (dump_frames) {
#>                 message("Loading Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             else {
#>                 message("Projecting Weights ", seq_start, " to ", 
#>                   seq_end, " of ", Nim)
#>             }
#>             pre_stack_weight_list = foreach(i = seq_start:seq_end, 
#>                 .packages = c("Rwcs", "Rfits"), export = c("image_list", 
#>                   "dump_dir", "weight_image")) %dopar% {
#>                 if (weight_image[i]) {
#>                   if (dump_frames) {
#>                     Rfits_read_image(paste0(dump_dir, "/weight_warp_", 
#>                       i, ".fits"))
#>                   }
#>                   else {
#>                     if (inherits(image_list[[i]], "Rfits_pointer")) {
#>                       temp_weight = Rfits_create_image(matrix(0L, 
#>                         dim(image_list[[i]])[1], dim(image_list[[i]])[2]), 
#>                         image_list[[i]]$keyvalues)
#>                     }
#>                     else {
#>                       temp_weight = image_list[[i]]
#>                     }
#>                     if (inherits(weight_list[[i]], "Rfits_pointer")) {
#>                       weight_list[[i]]$header = FALSE
#>                       temp_weight$imDat[] = weight_list[[i]][, 
#>                         ]
#>                     }
#>                     else if (inherits(weight_list[[i]], "Rfits_image")) {
#>                       temp_weight$imDat[] = weight_list[[i]]$imDat
#>                     }
#>                     else {
#>                       temp_weight$imDat[] = weight_list[[i]]
#>                     }
#>                     return(propaneWarp(image_in = temp_weight, 
#>                       keyvalues_out = keyvalues_out, dim_out = dim_out, 
#>                       doscale = FALSE, dotightcrop = TRUE, keepcrop = keepcrop, 
#>                       warpfield = pre_stack_image_list[[i - seq_start + 
#>                         1L]]$warpfield, cores = cores_warp, ...))
#>                   }
#>                 }
#>                 else {
#>                   return(weight_list[[i]])
#>                 }
#>             }
#>         }
#>         else {
#>             pre_stack_weight_list = weight_list
#>         }
#>         if (is.null(pre_stack_inVar_list)) {
#>             message("Stacking Images ", seq_start, " to ", seq_end, 
#>                 " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_cpp(post_image = post_stack_image, 
#>                   post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         else {
#>             message("Stacking Images and InVar ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mask_clip = mask_clip + temp_mask_clip
#>                 mode(temp_mask_clip) = "logical"
#>                 if (dump_frames) {
#>                   Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                     "/mask_warp_", i, ".fits"))
#>                 }
#>                 if (weight_image[i]) {
#>                   pre_weight = pre_stack_weight_list[[i]]$imDat
#>                 }
#>                 else {
#>                   pre_weight = weight_list[[i]]
#>                 }
#>                 .stack_image_inVar_cpp(post_image = post_stack_image, 
#>                   post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>                   pre_image = pre_stack_image_list[[i]]$imDat, 
#>                   pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>                   pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                     "YCUTLO")]), post_mask = temp_mask_clip)
#>             }
#>         }
#>         if (keep_extreme_pix) {
#>             message("Calculating Extreme Pixels ", seq_start, 
#>                 " to ", seq_end, " of ", Nim)
#>             temp_mask_clip = matrix(0L, dim_im[1], dim_im[2])
#>             for (i in 1:Nbatch_sub) {
#>                 i_stack = seq_start + i - 1L
#>                 if (i > 1) {
#>                   temp_mask_clip[] = 0L
#>                 }
#>                 temp_mask_clip[unique(c(post_stack_DT[cold == 
#>                   i_stack, ID], post_stack_DT[hot == i_stack, 
#>                   ID]))] = 1L
#>                 if (clip_dilate > 0) {
#>                   temp_mask_clip = .dilate_R(temp_mask_clip, 
#>                     size = clip_dilate)
#>                 }
#>                 mode(temp_mask_clip) = "logical"
#>                 xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>                 ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>                 new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                   post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                   post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                   ysub] == FALSE)
#>                 post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>                 post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>             }
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#605: message("Restacking without clipped cold/hot pixels")
#> Restacking without clipped cold/hot pixels
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#606: if (is.null(pre_stack_inVar_list)) {
#>     message("Stacking Images ", seq_start, " to ", seq_end, " of ", 
#>         Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_cpp(post_image = post_stack_image, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>             offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> } else {
#>     message("Stacking Images and InVar ", seq_start, " to ", 
#>         seq_end, " of ", Nim)
#>     for (i in 1:Nim) {
#>         if (i > 1) {
#>             temp_mask_clip[] = 0L
#>         }
#>         temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], 
#>             post_stack_DT[hot == i, ID]))] = 1L
#>         if (clip_dilate > 0) {
#>             temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>         }
#>         mask_clip = mask_clip + temp_mask_clip
#>         mode(temp_mask_clip) = "logical"
#>         if (dump_frames) {
#>             Rfits_write_image(temp_mask_clip, paste0(dump_dir, 
#>                 "/mask_warp_", i, ".fits"))
#>         }
#>         if (weight_image[i]) {
#>             pre_weight = pre_stack_weight_list[[i]]$imDat
#>         }
#>         else {
#>             pre_weight = weight_list[[i]]
#>         }
#>         .stack_image_inVar_cpp(post_image = post_stack_image, 
#>             post_inVar = post_stack_inVar, post_weight = post_stack_weight, 
#>             pre_image = pre_stack_image_list[[i]]$imDat, pre_inVar = pre_stack_inVar_list[[i]]$imDat, 
#>             pre_weight_sexp = pre_weight, offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>                 "YCUTLO")]), post_mask = temp_mask_clip)
#>         if (keep_extreme_pix) {
#>             xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>             ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>             new_cold = which(pre_stack_image_list[[i]]$imDat < 
#>                 post_stack_cold[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             new_hot = which(pre_stack_image_list[[i]]$imDat > 
#>                 post_stack_hot[xsub, ysub] & temp_mask_clip[xsub, 
#>                 ysub] == FALSE)
#>             post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>             post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>         }
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#658: message("Stacking Images and InVar ", seq_start, " to ", seq_end, 
#>     " of ", Nim)
#> Stacking Images and InVar 1 to 8 of 8
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#659: for (i in 1:Nim) {
#>     if (i > 1) {
#>         temp_mask_clip[] = 0L
#>     }
#>     temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>         i, ID]))] = 1L
#>     if (clip_dilate > 0) {
#>         temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#>     }
#>     mask_clip = mask_clip + temp_mask_clip
#>     mode(temp_mask_clip) = "logical"
#>     if (dump_frames) {
#>         Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>             i, ".fits"))
#>     }
#>     if (weight_image[i]) {
#>         pre_weight = pre_stack_weight_list[[i]]$imDat
#>     }
#>     else {
#>         pre_weight = weight_list[[i]]
#>     }
#>     .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>         post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>         pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>         offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>             "YCUTLO")]), post_mask = temp_mask_clip)
#>     if (keep_extreme_pix) {
#>         xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>         ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>         new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>             ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>         post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>         post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#>     }
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: if (i > 1) {
#>     temp_mask_clip[] = 0L
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#661: temp_mask_clip[] = 0L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#663: temp_mask_clip[unique(c(post_stack_DT[cold == i, ID], post_stack_DT[hot == 
#>     i, ID]))] = 1L
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#665: if (clip_dilate > 0) {
#>     temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#666: temp_mask_clip = .dilate_R(temp_mask_clip, size = clip_dilate)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#668: mask_clip = mask_clip + temp_mask_clip
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#669: mode(temp_mask_clip) = "logical"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#671: if (dump_frames) {
#>     Rfits_write_image(temp_mask_clip, paste0(dump_dir, "/mask_warp_", 
#>         i, ".fits"))
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#675: if (weight_image[i]) {
#>     pre_weight = pre_stack_weight_list[[i]]$imDat
#> } else {
#>     pre_weight = weight_list[[i]]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#678: pre_weight = weight_list[[i]]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#680: .stack_image_inVar_cpp(post_image = post_stack_image, post_inVar = post_stack_inVar, 
#>     post_weight = post_stack_weight, pre_image = pre_stack_image_list[[i]]$imDat, 
#>     pre_inVar = pre_stack_inVar_list[[i]]$imDat, pre_weight_sexp = pre_weight, 
#>     offset = unlist(pre_stack_image_list[[i]]$keyvalues[c("XCUTLO", 
#>         "YCUTLO")]), post_mask = temp_mask_clip)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#690: if (keep_extreme_pix) {
#>     xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#>     ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#>     new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>         ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#>     post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#>     post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#691: xsub = pre_stack_image_list[[i]]$keyvalues$XCUTLO:pre_stack_image_list[[i]]$keyvalues$XCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#692: ysub = pre_stack_image_list[[i]]$keyvalues$YCUTLO:pre_stack_image_list[[i]]$keyvalues$YCUTHI
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#694: new_cold = which(pre_stack_image_list[[i]]$imDat < post_stack_cold[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#695: new_hot = which(pre_stack_image_list[[i]]$imDat > post_stack_hot[xsub, 
#>     ysub] & temp_mask_clip[xsub, ysub] == FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#697: post_stack_cold[xsub, ysub][new_cold] = pre_stack_image_list[[i]]$imDat[new_cold]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#698: post_stack_hot[xsub, ysub][new_hot] = pre_stack_image_list[[i]]$imDat[new_hot]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#993: if (return_all == FALSE) {
#>     pre_stack_exp_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#997: weight_sel = (post_stack_weight != 0L)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#999: if (is.null(pre_stack_inVar_list)) {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_weight[weight_sel]
#> } else {
#>     post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#>     post_stack_inVar[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1002: post_stack_image[weight_sel] = post_stack_image[weight_sel]/post_stack_inVar[weight_sel]
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1003: post_stack_inVar[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1006: if (keep_extreme_pix) {
#>     post_stack_cold[!weight_sel] = NA
#>     post_stack_hot[!weight_sel] = NA
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1007: post_stack_cold[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1008: post_stack_hot[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1011: post_stack_image[!weight_sel] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1014: if (return_all == FALSE) {
#>     pre_stack_image_list = NULL
#>     pre_stack_inVar_list = NULL
#>     pre_stack_exp_list = NULL
#>     pre_stack_weight_list = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1021: keyvalues_out$EXTNAME = "image"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1022: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1023: keyvalues_out$R_VER = R.version$version.string
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1024: keyvalues_out$PANE_VER = as.character(packageVersion("ProPane"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1025: keyvalues_out$RWCS_VER = as.character(packageVersion("Rwcs"))
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1027: image_out = Rfits_create_image(post_stack_image, keyvalues = keyvalues_out, 
#>     keypass = FALSE, history = "Stacked with propaneStackWarpInVar")
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1032: keyvalues_out$EXTNAME = "weight"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1033: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1034: weight_out = Rfits_create_image(post_stack_weight, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1038: if (!is.null(post_stack_inVar)) {
#>     keyvalues_out$EXTNAME = "inVar"
#>     keyvalues_out$MAGZERO = magzero_out
#>     inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     inVar_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1039: keyvalues_out$EXTNAME = "inVar"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1040: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1041: inVar_out = Rfits_create_image(post_stack_inVar, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1048: if (!is.null(post_stack_exp)) {
#>     keyvalues_out$EXTNAME = "exp"
#>     keyvalues_out$MAGZERO = NULL
#>     exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     exp_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1049: keyvalues_out$EXTNAME = "exp"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1050: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1051: exp_out = Rfits_create_image(post_stack_exp, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1058: if (keep_extreme_pix) {
#>     post_stack_cold[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "cold"
#>     keyvalues_out$MAGZERO = magzero_out
#>     cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#>     post_stack_hot[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "hot"
#>     hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> } else {
#>     cold_out = NULL
#>     hot_out = NULL
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1059: post_stack_cold[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1061: keyvalues_out$EXTNAME = "cold"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1062: keyvalues_out$MAGZERO = magzero_out
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1063: cold_out = Rfits_create_image(post_stack_cold, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1067: post_stack_hot[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1069: keyvalues_out$EXTNAME = "hot"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1070: hot_out = Rfits_create_image(post_stack_hot, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1078: if (!is.null(mask_clip)) {
#>     mask_clip[weight_out$imDat == 0L] = NA
#>     keyvalues_out$EXTNAME = "clip"
#>     keyvalues_out$MAGZERO = NULL
#>     mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>         keypass = FALSE)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1079: mask_clip[weight_out$imDat == 0L] = NA
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1081: keyvalues_out$EXTNAME = "clip"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1082: keyvalues_out$MAGZERO = NULL
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1083: mask_clip = Rfits_create_image(mask_clip, keyvalues = keyvalues_out, 
#>     keypass = FALSE)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1088: time_taken = proc.time()[3] - timestart
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1089: message("Time taken: ", signif(time_taken, 4), " seconds")
#> Time taken: 16.79 seconds
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1091: if (return_all) {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>         exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> } else {
#>     output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>         exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>         which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>         dump_dir = dump_dir)
#> }
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1092: output = list(image = image_out, weight = weight_out, inVar = inVar_out, 
#>     exp = exp_out, cold = cold_out, hot = hot_out, clip = mask_clip, 
#>     image_pre_stack = pre_stack_image_list, inVar_pre_stack = pre_stack_inVar_list, 
#>     exp_pre_stack = pre_stack_exp_list, weight_pre_stack = pre_stack_weight_list, 
#>     which_overlap = which_overlap, time = time_taken, Nim = Nim, 
#>     dump_dir = dump_dir)
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1120: class(output) = "ProPane"
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1121: stopImplicitCluster()
#> debug at /Users/aaron/Rprojects/ProPane/R/propaneStackWarpInVar.R#1122: return(invisible(output))
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)

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 
#> -0.233271200628739 -0.030987709508136 -0.000203095969232  0.001236589012593 
#>            3rd Qu.               Max. 
#>  0.031376355397466  0.232873942686488
#> Pop Std Dev: 0.049925
#> MAD: 0.046216
#> Half 16-84 Quan (1s): 0.046615
#> Half 02-98 Quan (2s): 0.10092
#> Using 4386100 out of 4840000 (90.62%) data points (1305 < xlo & 48633 > 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 
#> -0.207729005007536 -0.027633839250850 -0.000317223281038  0.001115475877999 
#>            3rd Qu.               Max. 
#>  0.027797386074493  0.207101475624127
#> Pop Std Dev: 0.044753
#> MAD: 0.041077
#> Half 16-84 Quan (1s): 0.041484
#> Half 02-98 Quan (2s): 0.090567
#> Using 4378848 out of 4840000 (90.47%) data points (2314 < xlo & 54874 > 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 = cores, multitype = multitype,
                                 pattern = glob2rx('*image*fits') #to only stack images
                                 )
#> Time taken: 2.739 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.