Better recursive descent based Fibre Apertures. Correctly computes the fractional coverage down to a nominal recusion depth.

Author: asgr

DESCRIPTION

@@ -1,15 +1,15 @@
 Package: ProFound
 Type: Package
 Title: Photometry Tools
-Version: 1.26.3
-Date: 2025-10-09
+Version: 1.27.0
+Date: 2025-10-14
 Author: Aaron Robotham
 Maintainer: Aaron Robotham <aaron.robotham@uwa.edu.au>
 Description: Core package containing all the tools for simple and advanced source extraction. This is used to create inputs for 'ProFit', or for source detection, extraction and photometry in its own right.
 License: LGPL-3
 Depends: R (>= 3.1), Rfits (>= 1.8.0), magicaxis (>= 2.0.8), Rcpp (>= 1.0.2)
 Imports: data.table, celestial (>= 1.4.1), foreach, matrixStats, doParallel
-Suggests: ProFit, knitr, rmarkdown, EBImage, imager, LaplacesDemon, Rfast, Rfast2, fastmatch, snow, doSNOW, bigmemory, mvtnorm, Rwcs, Highlander (>= 0.1.7), ProPane
+Suggests: ProFit, knitr, rmarkdown, EBImage, imager, LaplacesDemon, Rfast, Rfast2, fastmatch, snow, doSNOW, bigmemory, mvtnorm, Rwcs, Highlander (>= 0.1.7), ProPane, plotrix
 VignetteBuilder: knitr
 LinkingTo: Rcpp
 NeedsCompilation: yes

R/RcppExports.R

@@ -16,6 +16,10 @@
     invisible(.Call(`_ProFound_interpolateLinearGrid`, xseq, yseq, tempmat_sky, output))
 }
 
+profoundAperCover <- function(x, y, cx, cy, radius, depth = 4L) {
+    .Call(`_ProFound_profoundAperCover`, x, y, cx, cy, radius, depth)
+}
+
 .dilate_cpp <- function(segim, kern, expand = 0L) {
     .Call(`_ProFound_dilate_cpp`, segim, kern, expand)
 }

R/profoundAperPhot.R

@@ -1,77 +1,84 @@
-.fluxcalcapp = function(x=NULL, y=NULL, rad2=NULL, flux=NULL, xcen=NA, ycen=NA, rad_app=NULL, centype='mean'){
-  if(is.null(rad2)){
-    if(is.na(xcen)){
-      if(centype == 'wt' | centype == 'mean'){
-        xcen = .meanwt(x, flux)
-      }else if(centype == 'max'){
-        xcen = x[which.max(flux)]
-      }else{
-        stop('centype must be max or mean!')
-      }
-    }
-    
-    if(is.na(ycen)){
-      if(centype == 'wt' | centype == 'mean'){
-        ycen = .meanwt(y, flux)
-      }else if(centype == 'max'){
-        ycen = y[which.max(flux)]
-      }else{
-        stop('centype must be max or mean!')
-      }
-    }
-    
-    if(xcen == 0 & ycen == 0){
-      rad2 = x^2 + y^2
+.fluxcalcapp = function(x=NULL, y=NULL, flux=NULL, xcen=NA, ycen=NA, rad_app=NULL, centype='mean', depth=4){
+  if(is.na(xcen)){
+    if(centype == 'wt' | centype == 'mean'){
+      xcen = .meanwt(x, flux)
+    }else if(centype == 'max'){
+      xcen = x[which.max(flux)]
     }else{
-      rad2 = (x - xcen)^2 + (y - ycen)^2 
+      stop('centype must be max or mean!')
     }
-  }else{
-    xcen = 0
-    ycen = 0
   }
 
-  if(is.na(rad_app)){
-    #left here for a different mode I tested for applying rad2 selection outside of this function.
-    #that wasn't any faster, so probably remove this if case in the future
-    Nsel = length(rad2)
-    rad_out = max(rad2)
-    sel_out = which(rad2 == rad_out)
-    
-    flux_app = sum(flux, na.rm=TRUE)
-    flux_min = mean(flux[sel_out], na.rm=TRUE)
-  }else{
-    sel = which(rad2 <= rad_app^2)
-    Nsel = length(sel)
-    
-    if(length(Nsel) == 0){
-      flux_app = NA_real_
-      flux_min = 0
+  if(is.na(ycen)){
+    if(centype == 'wt' | centype == 'mean'){
+      ycen = .meanwt(y, flux)
+    }else if(centype == 'max'){
+      ycen = y[which.max(flux)]
     }else{
-      if(all(is.na(flux[sel]))){
-        flux_app = NA_real_
-        flux_min = 0
-      }else{
-        rad_out = max(rad2[sel])
-        sel_out = which(rad2 == rad_out)
-        
-        flux_app = sum(flux[sel], na.rm=TRUE)
-        flux_min = mean(flux[sel_out], na.rm=TRUE)
-      }
+      stop('centype must be max or mean!')
     }
   }
 
+  # if(xcen == 0 & ycen == 0){
+  #   rad2 = x^2 + y^2
+  # }else{
+  #   rad2 = (x - xcen)^2 + (y - ycen)^2 
+  # }
+  
+  # if(is.na(rad_app)){
+  #   #left here for a different mode I tested for applying rad2 selection outside of this function.
+  #   #that wasn't any faster, so probably remove this if case in the future
+  #   Nsel = length(rad2)
+  #   rad_out = max(rad2)
+  #   sel_out = which(rad2 == rad_out)
+  #   
+  #   flux_app = sum(flux, na.rm=TRUE)
+  #   #flux_min = mean(flux[sel_out], na.rm=TRUE)
+  # }else{
+  #   sel = which(rad2 <= rad_app^2)
+  #   Nsel = length(sel)
+  #   
+  #   if(length(Nsel) == 0){
+  #     flux_app = NA_real_
+  #     #flux_min = 0
+  #   }else{
+  #     if(all(is.na(flux[sel]))){
+  #       flux_app = NA_real_
+  #       #flux_min = 0
+  #     }else{
+  #       rad_out = max(rad2[sel])
+  #       sel_out = which(rad2 == rad_out)
+  #       
+  #       flux_app = sum(flux[sel], na.rm=TRUE)
+  #       #flux_min = mean(flux[sel_out], na.rm=TRUE)
+  #     }
+  #   }
+  # }
+  
+  aper_frac = profoundAperCover(x, y, xcen, ycen, rad_app, depth=depth)
+  flux_app = sum(flux * aper_frac, na.rm=TRUE)
+  N_app = sum((!is.na(aper_frac))*aper_frac, na.rm=TRUE)
+  
+  suppressWarnings({
+    if(depth > 0){
+      flux_min = min(flux[aper_frac > 0 & aper_frac < 1 & flux > 0], na.rm=TRUE)
+    }else{
+      flux_min = min(flux[aper_frac > 0 & aper_frac <= 1 & flux > 0], na.rm=TRUE)
+    }
+  })
+  
   if(!isTRUE(is.finite(flux_min))){ #this catches NA, NaN, NULL, Inf events
     flux_min = 0 #don't want to penalise when masked or other weird events
   }else if(flux_min < 0){
     flux_min = 0 #don't want to penalise when in the sky noise 
   }
 
-  return(list(flux_app=flux_app, flux_min=flux_min, N=Nsel))
+  return(list(flux_app=flux_app, N_app=N_app, flux_min=flux_min))
 }
 
 profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyvalues=NULL, tar=NULL,
-                           pixscale=1, magzero=0, correction=TRUE, centype='mean', fluxtype='Raw',
-                           verbose=FALSE){
+                           pixscale=1, magzero=0, correction=TRUE, centype='mean',
+                           fluxtype='Raw', depth=4, verbose=FALSE){
   if(!is.null(image)){
     if(inherits(image, 'Rfits_image')){
       keyvalues = image$keyvalues
@@ -128,7 +135,7 @@ profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyva
     stop('fluxtype must be Jansky / Microjansky / Raw!')
   }
 
-  segID = x = y = flux = j = rad2 = NULL
+  segID = x = y = flux = j = NULL
 
   Rapp = (app_diam / 2 / pixscale) #in pixels
   Aapp = (pi * Rapp^2) #in pixels
@@ -207,7 +214,7 @@ profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyva
   match_segID = match(tempDT$segID, tar$segID)
   tempDT[, x:= x - tar[match_segID, 'xcen']]
   tempDT[, y:= y - tar[match_segID, 'ycen']]
-  tempDT[, rad2:= x^2 + y^2]
+  #tempDT[, rad2:= x^2 + y^2]
 
   #setkey(tempDT, segID, rad2)
 
@@ -218,7 +225,7 @@ profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyva
     if(verbose){message('  Aperture: ', app_diam[j], ' [asec]')}
     # the top one can completely remove segments in some weird edge cases, and doesn't appear to be faster. Use the second!
     #temp_app = tempDT[rad2 <= Rapp[j]^2, .fluxcalcapp(x=x, y=y, rad2=rad2, flux=flux, xcen=0, ycen=0, rad_app=NA), by=segID]
-    temp_app = tempDT[, .fluxcalcapp(rad2=rad2, flux=flux, xcen=0, ycen=0, rad_app=Rapp[j]), by=segID]
+    temp_app = tempDT[, .fluxcalcapp(x=x, y=y, flux=flux, xcen=0, ycen=0, rad_app=Rapp[j], depth=depth), by=segID]
 
     if(correction){
       temp_app$flux_app = temp_app$flux_app - (temp_app$N - Aapp[j])*temp_app$flux_min
@@ -228,9 +235,10 @@ profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyva
     return(data.frame(flux_app = temp_app$flux_app*fluxscale,
                       mag_app = mag_app,
                       SB_app = mag_app + 2.5*log10(Aapp[j]) + 5*log10(pixscale),
-                      N_app = temp_app$N,
-                      frac_app = temp_app$N/Aapp[j],
-                      flux_min = temp_app$flux_min*fluxscale)
+                      N_app = temp_app$N_app,
+                      frac_app = temp_app$N_app/Aapp[j],
+                      flux_min = temp_app$flux_min*fluxscale
+                      )
     )
   }
 
@@ -242,7 +250,8 @@ profoundAperPhot = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, keyva
 }
 
 profoundAperRan = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, Nran=100, keyvalues=NULL,
-                          pixscale=1, magzero=0, correction=TRUE, fluxtype='Raw', verbose=FALSE){
+                          pixscale=1, magzero=0, correction=TRUE, fluxtype='Raw', depth=4,
+                          verbose=FALSE){
   if(!is.null(image)){
     if(inherits(image, 'Rfits_image')){
       keyvalues = image$keyvalues
@@ -318,8 +327,8 @@ profoundAperRan = function(image=NULL, segim=NULL, app_diam=1, mask=NULL, Nran=1
   segim_ran[segim > 0L] = 0L
 
   output = profoundAperPhot(image=image, segim=segim_ran, app_diam=app_diam, keyvalues=keyvalues,
-                            tar=tar, pixscale=pixscale, magzero=magzero, correction=correction, 
-                            fluxtype=fluxtype, verbose=verbose)
+                            tar=tar, pixscale=pixscale, magzero=magzero, correction=correction,
+                            fluxtype=fluxtype, depth=depth, verbose=verbose)
 
   i = NULL
 

R/profoundSegim.R

@@ -95,26 +95,26 @@
     N100seg=length(flux)
   }
 
-  if(!is.na(Napp)){
-    Nsel = N100seg:(N100seg - Napp + 1)
-    Nsel = Nsel[Nsel>0]
-  }else{
-    Nsel = 0
-  }
+  # if(!is.na(Napp)){
+  #   Nsel = N100seg:(N100seg - Napp + 1)
+  #   Nsel = Nsel[Nsel>0]
+  # }else{
+  #   Nsel = 0
+  # }
 
   if(N100seg > 0){
 
     sumflux=sum(flux[good])
 
-    if(length(Nsel)>0){
-      if(Nsel[1]>0){
-        sumflux_app=sum(flux[good][Nsel])
-      }else{
-        sumflux_app=0
-      }
-    }else{
-      sumflux_app=0
-    }
+    # if(length(Nsel)>0){
+    #   if(Nsel[1]>0){
+    #     sumflux_app=sum(flux[good][Nsel])
+    #   }else{
+    #     sumflux_app=0
+    #   }
+    # }else{
+    #   sumflux_app=0
+    # }
 
     temp=cumsum(flux[good])/sumflux
 
@@ -139,21 +139,21 @@
 
   }else{
     sumflux=NA
-    sumflux_app=NA
+    #sumflux_app=NA
     N100seg=length(flux)
     N50seg=N100seg*0.5
     N90seg=N100seg*0.9
     cenfrac=NA
   }
 
   mode(sumflux)='numeric'
-  mode(sumflux_app)='numeric'
+  #mode(sumflux_app)='numeric'
   mode(N50seg)='numeric'
   mode(N90seg)='numeric'
   mode(cenfrac)='numeric'
   mode(cenfrac)='numeric'
 
-  return(list(flux=sumflux, flux_app=sumflux_app, N50seg=N50seg, N90seg=N90seg, N100seg=N100seg, cenfrac=cenfrac))
+  return(list(flux=sumflux, N50seg=N50seg, N90seg=N90seg, N100seg=N100seg, cenfrac=cenfrac))
 }
 
 .fluxcalcmin=function(flux){
@@ -791,15 +791,17 @@ profoundSegimStats=function(image=NULL, segim=NULL, mask=NULL, sky=NULL, skyRMS=
 
   x=NULL; y=NULL; flux=NULL; sky=NULL; skyRMS=NULL
 
-  fluxout = tempDT[, .fluxcalc(flux, Napp = NA), by = segID]
+  fluxout = tempDT[, .fluxcalc(flux), by = segID]
   #old very rough fibre mag stuff is being ignored (hence the NA above and commented out line below)
   #fluxout$flux_app[which(fluxout$flux_app > fluxout$flux)] = fluxout$flux[which(fluxout$flux_app > fluxout$flux)]
 
   if(!is.na(Rapp)){
     #newer more accurate fibre mag calculation
     temp_app = tempDT[, .fluxcalcapp(x=x, y=y, flux=flux, rad_app=Rapp), by=segID]
     #Here we correct by the lowest value pixel in the outer aperture
+    #no need now
     fluxout$flux_app = temp_app$flux_app - (temp_app$N - Aapp)*temp_app$flux_min
+    fluxout$flux_app = temp_app$flux_app
   }
 
   mag = profoundFlux2Mag(flux = fluxout$flux, magzero = magzero)

man/profoundAperCover.Rd

@@ -0,0 +1,73 @@
+\name{profoundAperCover}
+\alias{profoundAperCover}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{
+Fraction of Pixel in Aperture
+}
+\description{
+Calculates the fraction of pixels inside an aperture using recursive descent.
+}
+\usage{
+profoundAperCover(x, y, cx, cy, radius, depth = 4L)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{
+Numeric vector; the pixel x centre positions (where pixel centres are half integer in \code{ProFound}).
+}
+  \item{y}{
+Numeric vector; the pixel y centre positions (where pixel centres are half integer in \code{ProFound}).
+}
+  \item{cx}{
+Numeric scalar; the x centre of the aperture (in pixels).
+}
+  \item{cy}{
+Numeric scalar; the y centre of the aperture (in pixels).
+}
+  \item{radius}{
+Numeric scalar; the radius of the aperture (in pixels).
+}
+  \item{depth}{
+Integer scalar; the recursion depth to use.
+}
+}
+\details{
+Re the recursion \option{depth}, higher numbers are more accurate, where you probably do not want to set it below 3. Above 5 will suffer notable slow-down. 4 is usually about the performance versus accuracy sweet spot.
+}
+\value{
+Numeric vector the same length as \option{x} and \option{y}, which is what fraction of the pixel lies under the desired aperture.
+}
+\author{
+Aaron Robotham
+}
+
+\seealso{
+\code{\link{profoundAperPhot}}
+}
+\examples{
+library(magicaxis)
+library(plotrix)
+
+temp = expand.grid(1:10 - 0.5, 1:10 - 0.5)
+
+output = matrix(profoundAperCover(temp[,1], temp[,2], 5, 5, 3, depth=4), 10, 10)
+magimage(output, magmap=FALSE)
+draw.circle(5, 5, radius=3.2, border='red')
+
+output = matrix(profoundAperCover(temp[,1], temp[,2], 3.2, 3.2, 1.8, depth=4), 10, 10)
+magimage(output, magmap=FALSE)
+draw.circle(3.2, 3.2, radius=1.8, border='red')
+
+output = matrix(profoundAperCover(temp[,1], temp[,2], 7.5, 7.5, 5.8, depth=4), 10, 10)
+magimage(output, magmap=FALSE)
+draw.circle(7.5, 7.5, radius=5.8, border='red')
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory (show via RShowDoc("KEYWORDS")):
+% \keyword{ ~kwd1 }
+% \keyword{ ~kwd2 }
+% Use only one keyword per line.
+% For non-standard keywords, use \concept instead of \keyword:
+% \concept{ ~cpt1 }
+% \concept{ ~cpt2 }
+% Use only one concept per line.