R/profoundProFound.R
In asgr/ProFound: Photometry Tools
Defines functions
.selectCoG
plot.profound
profoundProFound
Documented in
plot.profound
profoundProFound
profoundProFound=function(image=NULL, segim=NULL, objects=NULL, mask=NULL, skycut=1.5, pixcut=3,
tolerance=4, ext=2, reltol=0, cliptol=Inf, sigma=1, smooth=TRUE,
SBlim=NULL, SBdilate=2, SBN100=100, size=5, shape='disc', iters=6,
threshold=1.05, magzero=0, gain=NULL, pixscale=1, sky=NULL, skyRMS=NULL,
redosegim=FALSE, redosky=TRUE, redoskysize=21, box=c(100,100), grid=box,
skygrid_type = 'new', type='bicubic', skytype='median', skyRMStype='quanlo', roughpedestal=FALSE,
sigmasel=1, skypixmin=prod(box)/2, boxadd=box/2, boxiters=0, conviters=100, iterskyloc=FALSE,
deblend=FALSE, df=3, radtrunc=2, iterative=FALSE, doChiSq=FALSE, doclip=TRUE,
shiftloc = FALSE, paddim = TRUE, keyvalues=NULL, verbose=FALSE,
plot=FALSE, stats=TRUE, rotstats=FALSE, boundstats=FALSE,
nearstats=boundstats, groupstats=boundstats, group=NULL,
groupby='segim_orig', offset=1, haralickstats=FALSE, sortcol="segID",
decreasing=FALSE, lowmemory=FALSE, keepim=TRUE, watershed='ProFound',
pixelcov=FALSE, deblendtype='fit', psf=NULL, fluxweight='sum',
convtype = 'brute', convmode = 'extended', fluxtype='Raw',
app_diam=1, Ndeblendlim=Inf, static_photom = FALSE, rem_mask=FALSE, ...){
if(verbose){message('Running ProFound:')}
timestart=proc.time()[3]
call=match.call(expand.dots=TRUE)
if(length(box)==1){
box=rep(box,2)
if(missing(grid)){grid=box}
if(missing(boxadd)){boxadd=box/2}
if(missing(skypixmin)){skypixmin=prod(box)/2}
}
if(length(grid)==1){
grid=rep(grid,2)
}
if(length(boxadd)==1){
boxadd=rep(boxadd,2)
}
fluxtype = tolower(fluxtype)
if(fluxtype=='raw' | fluxtype=='adu' | fluxtype=='adus'){
if(verbose){message('Using raw flux units')}
fluxscale=1
}else if (fluxtype=='jansky'){
if(verbose){message('Using Jansky flux units (WARNING: magzero must take system to AB)')}
fluxscale=10^(-0.4*(magzero-8.9))
}else if (fluxtype=='microjansky'){
if(verbose){message('Using Micro-Jansky flux units (WARNING: magzero must take system to AB)')}
fluxscale=10^(-0.4*(magzero-23.9))
}else{
stop('fluxtype must be Jansky / Microjansky / Raw!')
}
#Split out image and keyvalues parts of input:
if(!is.null(image)){
if(inherits(image, 'Rfits_image')){
keyvalues = image$keyvalues
image = image$imDat
}else if(inherits(image, 'Rfits_pointer')){
keyvalues = image$keyvalues
image = image[,]$imDat
}else if(inherits(image, 'matrix')){
'Do nothing'
}else{
stop('As of ProFound v1.21.0 only Rfits_image FITS inputs are allowed. Please install from GitHub asgr/Rfits')
}
}else{
stop('Missing image - this is a required input!')
}
if(!is.null(keyvalues)){
if(!inherits(keyvalues, 'Rfits_keylist')){
if(is.list(keyvalues)){
class(keyvalues) = 'Rfits_keylist'
}else{
stop('keyvalues is the wrong format- should be a list!')
}
}
}
if(box[1] > ceiling(dim(image)[1]/3)){
box[1] = ceiling(dim(image)[1]/3)
message('dim(image)[1]/box[1] must be >=3, box[1] modified to ',box[1])
}
if(box[2] > ceiling(dim(image)[2]/3)){
box[2] = ceiling(dim(image)[2]/3)
message('dim(image)[2]/box[2] must be >=3, box[2] modified to ',box[2])
}
if(grid[1] > ceiling(dim(image)[1]/3)){
grid[1] = ceiling(dim(image)[1]/3)
message('dim(image)[1]/grid[1] must be >=3, grid[1] modified to ',grid[1])
}
if(grid[2] > ceiling(dim(image)[2]/3)){
grid[2] = ceiling(dim(image)[2]/3)
message('dim(image)[2]/grid[2] must be >=3, grid[2] modified to ',grid[2])
}
if(verbose){message(paste('Supplied image is',dim(image)[1],'x',dim(image)[2],'pixels'))}
#Treat image NAs as masked regions:
badpix=NULL
if(!is.null(mask)){
mask=mask*1L #Looks silly, but this ensures a logical mask becomes integer.
if(length(mask)==1 & !is.na(mask[1])){
maskflag=mask
mask=matrix(0L,dim(image)[1],dim(image)[2])
mask[image==maskflag]=1L
}
if(anyNA(image)){
badpix=which(is.na(image))
mask[badpix]=1L
image[badpix]=0
}
if(is.numeric(mask)){
mask=as.integer(mask)
attributes(mask)$dim = dim(image)
}
}else{
if(anyNA(image)){
mask=matrix(0L,dim(image)[1],dim(image)[2])
badpix=which(is.na(image))
mask[badpix]=1L
image[badpix]=0
}
}
#if(!is.null(segim) & !is.null(mask)){
# segim=segim*(1-mask) #I don't think we actually need this
#}
#Get the pixel scale, if possible and not provided:
if(missing(pixscale) & !is.null(keyvalues)){
pixscale = pixscale(keyvalues)
if(verbose){message(paste('Extracted pixel scale from keyvalues provided:',signif(pixscale,4),'asec/pixel'))}
}else{
if(verbose){message(paste('Using suggested pixel scale:',signif(pixscale,4),'asec/pixel'))}
}
imarea = prod(dim(image))*pixscale^2/(3600^2)
if(verbose){message(paste('Supplied image is',signif(dim(image)[1]*pixscale/60,4),'x',signif(dim(image)[2]*pixscale/60,4),'amin, ', signif(imarea,4),'deg-sq'))}
if(is.null(objects)){
if(!is.null(segim)){
objects = matrix(0L,dim(segim)[1],dim(segim)[2])
objects[] = as.logical(segim)
}
}#else{
# objects = objects*1L #Looks silly, but this ensures a logical mask becomes integer.
# }
if(static_photom & !is.null(segim)){
#set all the flags to make static photom
sky = 0
iters = 0
redosky = FALSE
redosegim = FALSE
}
#Check for user provided sky, and compute if missing:
hassky = !is.null(sky)
hasskyRMS = !is.null(skyRMS)
if((hassky==FALSE | hasskyRMS==FALSE) & is.null(segim)){
if(hassky==FALSE){
sky = 0
}
if(verbose){message(paste('Making initial sky map -',round(proc.time()[3]-timestart,3),'sec'))}
roughsky = profoundMakeSkyGrid(image=image, objects=objects, mask=mask, sky=sky, box=box,
grid=grid, skygrid_type=skygrid_type, type=type,
skytype=skytype, skyRMStype=skyRMStype, sigmasel=sigmasel,
skypixmin=skypixmin, boxadd=boxadd, boxiters=0, conviters=conviters,
doclip=doclip, shiftloc=shiftloc, paddim=paddim, rem_mask=rem_mask)
if(roughpedestal){
roughsky$sky = median(roughsky$sky,na.rm=TRUE)
roughsky$skyRMS = median(roughsky$skyRMS,na.rm=TRUE)
}
if(hassky==FALSE){
sky = roughsky$sky
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
}
if(hasskyRMS==FALSE){
skyRMS = roughsky$skyRMS
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}
rm(roughsky)
}else{
if(verbose){message("Skipping making initial sky map - User provided sky and sky RMS, or user provided segim")}
}
#Make the initial segmentation map, if not provided.
if(is.null(segim)){
if(verbose){message(paste('Making initial segmentation image -',round(proc.time()[3]-timestart,3),'sec'))}
segim = profoundMakeSegim(image=image, objects=objects, mask=mask, tolerance=tolerance,
ext=ext, reltol=reltol, cliptol=cliptol, sigma=sigma,
smooth=smooth, pixcut=pixcut, skycut=skycut, SBlim=SBlim,
sky=sky, skyRMS=skyRMS, magzero=magzero, pixscale=pixscale,
verbose=verbose, watershed=watershed, plot=FALSE, stats=FALSE)
objects = segim$objects
segim = segim$segim
}else{
redosegim = FALSE
if(verbose){message("Skipping making an initial segmentation image - User provided segim")}
}
if(any(segim>0)){
if((hassky==FALSE | hasskyRMS==FALSE)){
if(redosky){
if(verbose){message(paste('Doing initial aggressive dilation -',round(proc.time()[3]-timestart,3),'sec'))}
objects_redo=profoundMakeSegimDilate(segim=objects, size=redoskysize, shape=shape, sky=sky, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$objects
}else{
objects_redo=objects
}
if(verbose){message(paste('Making better sky map -',round(proc.time()[3]-timestart,3),'sec'))}
if(hasskyRMS==FALSE){rm(skyRMS)}
bettersky = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky,
box=box, skygrid_type=skygrid_type, grid=grid,
type=type, skytype=skytype, skyRMStype=skyRMStype,
sigmasel=sigmasel, skypixmin=skypixmin, boxadd=boxadd,
boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)
if(hassky==FALSE){
sky = bettersky$sky
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
}
if(hasskyRMS==FALSE){
skyRMS = bettersky$skyRMS
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}
rm(bettersky)
if(redosegim){
if(verbose){message(paste('Making better segmentation image -',round(proc.time()[3]-timestart,3),'sec'))}
imagescale = (image-sky)/skyRMS
imagescale[!is.finite(imagescale)] = 0
if(!is.null(SBlim) & !missing(magzero)){
imagescale[imagescale<skycut | sky<profoundSB2Flux(SBlim, magzero, pixscale)]=0
}else{
imagescale[imagescale<skycut] = 0
}
if(!is.null(mask)){
imagescale[mask!=0] = 0
}
segim[imagescale==0] = 0
objects[segim==0] = 0
rm(imagescale)
}
}else{
if(verbose){message("Skipping making better sky map - User provided sky and sky RMS")}
}
if(iters>0 | iterskyloc){
if(verbose){message(paste('Calculating initial segstats -',round(proc.time()[3]-timestart,3),'sec'))}
# if(length(sky)>1){
# skystats = .profoundFluxCalcMin(image = sky - median(sky,na.rm=TRUE), segim=segim, mask=mask) #run on sky
# skystats = skystats$flux / skystats$N100 #get per pixel mean flux per segment
# #skymed = median(skystats, na.rm=TRUE) #median per pixel mean flux per segment for the whole image
# }else{
# skystats = 0 #specified
# #skymed = sky #specified
# }
image_sky = image - sky
segstats = .profoundFluxCalcMin(image=image_sky, segim=segim, mask=mask) #run on initial image
#segstats$flux = segstats$flux - (skystats * segstats$N100) #remove the local sky component
origfrac = segstats$flux #set to initial fluxes
segim_orig = segim
expand_segID = segstats[,'segID']
SBlast = rep(Inf,length(expand_segID))
selseg = rep(0,length(expand_segID))
if(is.null(SBdilate)){
SBdilate = Inf
}else{
skyRMSstats = .profoundFluxCalcMin(image=skyRMS, segim=segim, mask=mask) #run on sky
skyRMSstats = skyRMSstats$flux/skyRMSstats$N100 #get per pixel mean flux per segment
SBdilate = skyRMSstats * 10^(-0.4*SBdilate)
SBdilate[!is.finite(SBdilate)]=Inf
}
if(iters>0){
if(verbose){message('Doing dilations:')}
for(i in 1:(iters)){
if(verbose){message(paste('Iteration',i,'of',iters,'-',round(proc.time()[3]-timestart,3),'sec'))}
segim_new = profoundDilate(segim=segim, size=size, shape=shape, expand=expand_segID, iters=1)
#segim_new = profoundMakeSegimDilate(segim=segim, expand=expand_segID, size=size, shape=shape, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$segim #dilate
segstats_new = .profoundFluxCalcMin(image=image_sky, segim=segim_new, mask=mask) #run on image with dilated segments
N100diff = (segstats_new$N100 - segstats$N100)
SBnew = (segstats_new$flux - segstats$flux) / N100diff #calculate the surface brightness of the new grown anulus only
fluxgrowthratio = segstats_new$flux / segstats$flux #calculate flux growth ratio
#skyfrac = abs((skystats * N100diff) / (segstats_new$flux - segstats$flux)) #estimate how much of the flux growth might be coming from unusual local sky
expand_segID = segstats[which((fluxgrowthratio > threshold | (SBnew > SBdilate & N100diff>SBN100)) & segstats_new$flux>0 & SBnew < SBlast & selseg==(i-1)),'segID']
expand_segID = expand_segID[is.finite(expand_segID)] #safety first
if(length(expand_segID)==0){break}
updateID = which(segstats$segID %in% expand_segID)
selseg[updateID] = i #iteration number for segments flagged for growth
segstats[updateID,] = segstats_new[updateID,] #update only the segstats for things that are growing, the rest can be ignored
SBlast = SBnew
if('fastmatch' %in% .packages()){ #dilate segments that pass tests
selpix = which(fastmatch::fmatch(segim_new, expand_segID, nomatch = 0L) > 0)
}else{
selpix = which(segim_new %in% expand_segID)
}
segim[selpix] = segim_new[selpix]
}
}
origfrac = origfrac / segstats$flux
if(iterskyloc){
segim_skyloc = profoundMakeSegimDilate(segim=segim, size=size, shape=shape, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$segim
segstats_sky = .profoundFluxCalcMin(image=image, segim=segim_skyloc-segim, mask=mask)
#segstats$flux = segstats$flux + (skystats * segstats$N100) #add back the local sky component
skyseg_temp = segstats_sky$flux / segstats_sky$N100
skyseg_temp[!is.finite(skyseg_temp)] = 0
skyseg_mean = rep(NA,length(segstats$segID))
skyseg_mean[match(segstats_sky$segID,segstats$segID)] = skyseg_temp
}else{
skyseg_mean = NA
}
objects = matrix(0L,dim(segim)[1],dim(segim)[2])
objects[] = as.logical(segim)
}else{
if(verbose){message('Iters set to 0 - keeping segim un-dilated')}
segim_orig = segim
selseg = 0
origfrac = 1
skyseg_mean = NA
}
objects_redo = profoundMakeSegimDilate(segim=objects, mask=mask, size=redoskysize, shape=shape, sky=sky, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$objects
if(redosky){
if(redoskysize %% 2 == 0){redoskysize=redoskysize+1}
if(verbose){message(paste('Doing final aggressive dilation -',round(proc.time()[3]-timestart,3),'sec'))}
if(verbose){message(paste('Making final sky map -',round(proc.time()[3]-timestart,3),'sec'))}
rm(skyRMS)
sky = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type=skygrid_type, grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)$sky
if(verbose){message(paste('Making final sky RMS map -',round(proc.time()[3]-timestart,3),'sec'))}
skyRMS = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type=skygrid_type, grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)$skyRMS
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}else{
if(verbose){message("Skipping making final sky map - redosky set to FALSE")}
}
if(doChiSq){
if(verbose){message(paste('Making sky LL ChiSq map -',round(proc.time()[3]-timestart,3),'sec'))}
skyChiSqMap = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type='old', grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doChiSq=TRUE,
doclip=doclip, shiftloc=shiftloc, paddim=paddim, rem_mask=rem_mask)$skyChiSq
}else{
skyChiSqMap = NA
}
Norig = tabulate(segim_orig)
if(is.function(pixelcov)){
cor_err_func = pixelcov
}else{
if(pixelcov){
if(verbose){message(paste('Calculating pixel covariance -',round(proc.time()[3]-timestart,3),'sec'))}
if(is.null(profoundMakeSkyGrid)){
cor_err_func = profoundPixelCorrelation(image=image, objects=objects, mask=mask, sky=sky, skyRMS=skyRMS,
fft=FALSE, lag=apply(expand.grid(c(1,2,4),c(1,10,100,1000,1e4)),MARGIN=1,FUN=prod))$cor_err_func
}else{
cor_err_func = profoundPixelCorrelation(image=image, objects=objects_redo, mask=mask, sky=sky, skyRMS=skyRMS,
fft=FALSE, lag=apply(expand.grid(c(1,2,4),c(1,10,100,1000,1e4)),MARGIN=1,FUN=prod))$cor_err_func
}
}else{
cor_err_func = NULL
}
}
if(lowmemory){
image = image - sky
sky = 0
skyRMS = 0
segim_orig = NULL
objects = NULL
objects_redo = NULL
}
if(stats & !is.null(image)){
if(verbose){message(paste('Calculating final segstats for',length(which(tabulate(segim)>0)),'objects -',round(proc.time()[3]-timestart,3),'sec'))}
if(verbose){message(paste(' - magzero =', signif(magzero,4)))}
if(verbose){
if(is.null(gain)){
message(paste(' - gain = NULL (ignored)'))
}else{
message(paste(' - gain =', signif(gain,4)))
}
}
if(verbose){message(paste(' - pixscale =', signif(pixscale,4)))}
if(verbose){message(paste(' - rotstats =', rotstats))}
if(verbose){message(paste(' - boundstats =', boundstats))}
segstats = profoundSegimStats(image=image, segim=segim, mask=mask, sky=sky, skyRMS=skyRMS,
magzero=magzero, gain=gain, pixscale=pixscale, keyvalues=keyvalues,
sortcol=sortcol, decreasing=decreasing, rotstats=rotstats,
boundstats=boundstats, offset=offset, cor_err_func=cor_err_func,
app_diam=app_diam)
segstats = cbind(segstats, iter=selseg, origfrac=origfrac, Norig=Norig[segstats$segID], skyseg_mean=skyseg_mean)
segstats = cbind(segstats, flag_keep=segstats$origfrac >= median(segstats$origfrac[segstats$iter==iters]) | segstats$iter<iters)
}else{
if(verbose){message("Skipping segmentation statistics - segstats set to FALSE")}
segstats = NULL
}
if(nearstats){
near = profoundSegimNear(segim=segim, offset=offset)
}else{
near = NULL
}
if(deblend){
groupstats = TRUE
}
if(groupstats){
if(verbose){message(paste(' - groupstats = TRUE - ',round(proc.time()[3]-timestart,3),'sec'))}
if(groupby=='segim'){
if(is.null(group)){
group=profoundSegimGroup(segim)
}
}else if(groupby=='segim_orig'){
if(is.null(group)){
#message(round(proc.time()[3]-timestart,3))
group = profoundSegimGroup(segim_orig)
if(any(group$groupsegID$Ngroup>1)){
#message(round(proc.time()[3]-timestart,3))
group$groupim = profoundSegimKeep(segim=segim, segID_merge=group$groupsegID[group$groupsegID$Ngroup>1,'segID'])
#message(round(proc.time()[3]-timestart,3))
group$groupsegID$Npix = tabulate(group$groupim)[group$groupsegID$groupID]
}
}
}else{
stop('Non legal groupby option, must be segim or segim_orig!')
}
if(stats & !is.null(image) & !is.null(group)){
groupstats = profoundSegimStats(image=image, segim=group$groupim, mask=mask,
sky=sky, skyRMS=skyRMS, magzero=magzero, gain=gain,
pixscale=pixscale, keyvalues=keyvalues, sortcol=sortcol,
decreasing=decreasing, rotstats=rotstats, boundstats=boundstats,
offset=offset, cor_err_func=cor_err_func, app_diam=app_diam)
colnames(groupstats)[1] = 'groupID'
}else{
groupstats = NULL
}
}else{
if(verbose){message(' - groupstats = FALSE')}
group = NULL
groupstats = NULL
}
if(deblend & stats & !is.null(image) & any(group$groupsegID$Ngroup>1)){
if(verbose){message(paste(' - deblend = TRUE - ',round(proc.time()[3]-timestart,3),'sec'))}
tempblend = profoundFluxDeblend(image=image-sky, segim=segim, segstats=segstats,
groupim=group$groupim, groupsegID=group$groupsegID,
magzero=magzero, df=df, radtrunc=radtrunc, iterative=iterative,
doallstats=TRUE, deblendtype=deblendtype, psf=psf,
fluxweight=fluxweight, convtype=convtype, convmode=convmode,
Ndeblendlim = Ndeblendlim)
if(!is.null(tempblend)){
segstats = cbind(segstats,tempblend[,-2])
}
}else{
if(verbose){message(' - deblend = FALSE')}
}
if(haralickstats){
if(requireNamespace("EBImage", quietly = TRUE)){
scale = 10^(0.4*(30-magzero))
temphara = (image-sky)*scale
if(!is.null(mask)){
temphara[mask!=0] = 0
}
temphara[!is.finite(temphara)] = 0
haralick = as.data.frame(EBImage::computeFeatures.haralick(segim,temphara))
haralick = haralick[segstats$segID,]
}else{
if(verbose){
message('The EBImage package is needed to compute Haralick statistics.')
haralick = NULL
}
}
}else{
haralick = NULL
}
if(plot){
if(verbose){message(paste('Plotting segments -',round(proc.time()[3]-timestart,3),'sec'))}
if(any(is.finite(sky))){
profoundSegimPlot(image=image-sky, segim=segim, mask=mask, keyvalues=keyvalues, ...)
}else{
profoundSegimPlot(image=image, segim=segim, mask=mask, keyvalues=keyvalues, ...)
}
}else{
if(verbose){message("Skipping segmentation plot - plot = FALSE")}
}
if(is.null(SBlim)){
SBlim = NULL
}else if(is.numeric(SBlim)){
SBlimtemp = profoundFlux2SB(flux=skyRMS*skycut, magzero=magzero, pixscale=pixscale)
SBlimtemp = matrix(SBlimtemp,dim(skyRMS)[1],dim(skyRMS)[2])
SBlimtemp[which(SBlimtemp>SBlim)] = SBlim
SBlim = SBlimtemp
}else if(SBlim[1]=='get' & skycut> -Inf){
SBlim = profoundFlux2SB(flux=skyRMS*skycut, magzero=magzero, pixscale=pixscale)
}
if(is.null(keyvalues)){keyvalues = NULL}
if(keepim==FALSE){image = NULL; mask = NULL}
if(is.null(mask)){mask = NULL}
if(!is.null(badpix)){image[badpix] = NA}
row.names(segstats) = NULL
segstats[,grep('flux',colnames(segstats))] = fluxscale*segstats[,grep('flux',colnames(segstats))]
if(!is.null(groupstats)){
groupstats[,grep('flux',colnames(groupstats))] = fluxscale*groupstats[,grep('flux',colnames(groupstats))]
}
if(stats){
cutsky = (image - sky) / skyRMS
cutsky[!is.finite(cutsky)] = NA
if(!is.null(mask)){
cutsky[mask>0] = NA
}
if(!is.null(objects_redo)){
cutsky[objects_redo == 1] = NA
}else{
cutsky[objects == 1] = NA
}
cutsky = cutsky[which(cutsky<=0)] #Only consider sky pixels below the mode
if(length(cutsky) > 0){
df = length(cutsky) - 1
skyChiSq = sum(cutsky^2, na.rm =TRUE)
skyLL = dchisq(skyChiSq, df=df, log=TRUE)
skyChiSq = skyChiSq / df
}else{
skyLL = NA
skyChiSq = NA
}
}else{
skyLL = NULL
skyChiSq = NULL
}
if(verbose){message(paste('ProFound is finished! -',round(proc.time()[3]-timestart,3),'sec'))}
output = list(segim=segim, segim_orig=segim_orig, objects=objects, objects_redo=objects_redo,
sky=sky, skyRMS=skyRMS, image=image, mask=mask, segstats=segstats,
Nseg=dim(segstats)[1], near=near, group=group, groupstats=groupstats,
haralick=haralick, keyvalues=keyvalues, SBlim=SBlim, magzero=magzero, dim=dim(segim),
pixscale=pixscale, imarea=imarea, skyLL=skyLL, skyChiSq=skyChiSq, skyChiSqMap=skyChiSqMap,
gain=gain, call=call, date=date(), time=proc.time()[3]-timestart, ProFound.version=packageVersion('ProFound'),
R.version=R.version)
}else{
if(is.null(keyvalues)){keyvalues = NULL}
if(keepim==FALSE){image = NULL; mask = NULL}
if(is.null(mask)){mask = NULL}
if(!is.null(badpix)){image[badpix] = NA}
if(verbose){message('No objects in segmentation map - skipping dilations and CoG')}
if(verbose){message(paste('ProFound is finished! -',round(proc.time()[3]-timestart,3),'sec'))}
output = list(segim=NULL, segim_orig=NULL, objects=NULL, objects_redo=NULL, sky=sky,
skyRMS=skyRMS, image=image, mask=mask, segstats=NULL, Nseg=0, near=NULL,
group=NULL, groupstats=NULL, haralick=NULL, keyvalues=keyvalues, SBlim=NULL,
magzero=magzero, dim=dim(segim), pixscale=pixscale, imarea=imarea, skyLL=NULL, skyChiSq=NULL,
skyChiSqMap=NULL, gain=gain, call=call, date=date(), time=proc.time()[3]-timestart,
ProFound.version=packageVersion('ProFound'), R.version=R.version)
}
class(output) = 'profound'
return(invisible(output))
}
plot.profound=function(x, logR50=TRUE, dmag=0.5, hist='sky', ...){
suppressWarnings({
if(!inherits(x, 'profound')){
stop('Object class is not of type profound!')
}
if(is.null(x$image)){
stop('Missing image!')
}
if(is.null(x$segim)){
stop('Missing segmentation map!')
}
if(is.null(x$sky)){
x$sky=matrix(0, x$dim[1], x$dim[2])
}
if(length(x$sky)==1){
x$sky=matrix(x$sky, x$dim[1], x$dim[2])
}
if(is.null(x$skyRMS)){
x$skyRMS=matrix(1, x$dim[1], x$dim[2])
}
if(length(x$skyRMS)==1){
x$skyRMS=matrix(x$skyRMS, x$dim[1], x$dim[2])
}
segdiff = x$segim - x$segim_orig
segdiff[segdiff<0]=0
if(!is.null(x$mask)){
x$image[x$mask>0] = NA
x$mask[is.na(x$image)] = 1
}
image = (x$image-x$sky)/x$skyRMS
image[x$skyRMS < 0] = NA
if(is.null(x$objects)){
x$objects = matrix(0L,dim(x$segim)[1],dim(x$segim)[2])
x$objects[] = as.logical(x$segim)
}
#Masked and unmasked stats:
if(!is.null(x$mask)){
stat_mean_sky = signif(mean(x$sky[x$mask==0], na.rm=TRUE),4)
stat_sd_sky = signif(sd(x$sky[x$mask==0], na.rm=TRUE),4)
stat_mean_skyRMS = signif(mean(x$skyRMS[x$mask==0], na.rm=TRUE),4)
stat_sd_skyRMS = signif(sd(x$skyRMS[x$mask==0], na.rm=TRUE),4)
stat_cor_sky_skyRMS = signif(cor(as.numeric(x$sky[x$mask==0]), as.numeric(x$skyRMS[x$mask==0])^2, use="pairwise.complete.obs"),4)
if(!is.null(x$objects_redo)){
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$mask==0 & x$objects_redo>0]), as.numeric(x$image[x$mask==0 & x$objects_redo>0] - as.numeric(x$sky[x$mask==0 & x$objects_redo>0])), use="pairwise.complete.obs"),4)
}else{
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$mask==0 & x$objects>0]), as.numeric(x$image[x$mask==0 & x$objects>0] - as.numeric(x$sky[x$mask==0 & x$objects>0])), use="pairwise.complete.obs"),4)
}
}else{
stat_mean_sky = signif(mean(x$sky, na.rm=TRUE),4)
stat_sd_sky = signif(sd(x$sky, na.rm=TRUE),4)
stat_mean_skyRMS = signif(mean(x$skyRMS, na.rm=TRUE),4)
stat_sd_skyRMS = signif(sd(x$skyRMS, na.rm=TRUE),4)
stat_cor_sky_skyRMS = signif(cor(as.numeric(x$sky), as.numeric(x$skyRMS)^2, use="pairwise.complete.obs"),4)
if(!is.null(x$objects_redo)){
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$objects_redo>0]), as.numeric(x$image[x$objects_redo>0] - as.numeric(x$sky[x$objects_redo>0])), use="pairwise.complete.obs"),4)
}else{
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$objects>0]), as.numeric(x$image[x$objects>0] - as.numeric(x$sky[x$objects>0])), use="pairwise.complete.obs"),4)
}
}
if(requireNamespace("imager", quietly = TRUE)){
outline = x$objects - as.matrix(imager::erode(imager::as.cimg(x$objects),imager::imfill(3,3,val=1)))
}else{
outline = NULL
}
cmap = hcl.colors(100, 'RdYlBu', rev=TRUE)
maximg = quantile(abs(image[is.finite(image)]), 0.995, na.rm=TRUE)
stretchscale = 1/median(abs(image), na.rm=TRUE)
layout(matrix(1:9, 3, byrow=TRUE))
if(!is.null(x$keyvalues)){
par(mar=c(3.5,3.5,0.5,0.5)) #plot 1 image
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=image, keyvalues=x$keyvalues, stretchscale=stretchscale, locut=-maximg, hicut=maximg, range=c(-1,1), type='num', zlim=c(-1,1), col=cmap)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 2 seg
profoundSegimPlot(profound=x)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 3 dilate
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=image, keyvalues=x$keyvalues)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
magimage(segdiff, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE, add=TRUE)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 4 mag hist
if(is.null(x$imarea)){
imarea=prod(x$dim)*x$pixscale^2/(3600^2)
}else{
imarea=x$imarea
}
temphist=maghist(x$segstats$mag, log='y', scale=(2*dmag)/x$imarea, breaks=seq(floor(min(x$segstats$mag, na.rm = TRUE)), ceiling(max(x$segstats$mag, na.rm = TRUE)),by=0.5),
xlab='mag', ylab=paste('#/deg-sq/d',dmag,'mag',sep=''), grid=TRUE, verbose=FALSE)
#magplot(temphist, log='y', xlab='mag', ylab=expression('#'/'deg-sq'/'dmag'), grid=TRUE)
ymax=log10(max(temphist$counts,na.rm = T))
xmax=temphist$mids[which.max(temphist$counts)]
abline(ymax - xmax*0.4, 0.4, col='red')
abline(v=xmax+0.25, col='red')
axis(side=1, at=xmax+0.25, labels=xmax+0.25, tick=FALSE, line=-1, col.axis='red')
legend('topleft', legend=paste('N:',length(x$segstats$mag)), bty='n', text.col='black')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 5 sky
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=x$sky-median(x$sky,na.rm=TRUE), keyvalues=x$keyvalues, qdiff=TRUE)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky',paste0('Mean: ',stat_mean_sky),paste0('SD: ',stat_sd_sky)), bty='n', text.col='magenta')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 6 skyRMS
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=x$skyRMS, keyvalues=x$keyvalues)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'orange'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky RMS',paste0('Mean: ',stat_mean_skyRMS),paste0('SD: ',stat_sd_skyRMS)), bty='n', text.col='magenta')
#plot 7 sky den
if(hist=='iters'){
maghist(x$segstats$iter, breaks=seq(-0.5,max(x$segstats$iter, na.rm=TRUE)+0.5,by=1), majorn=max(x$segstats$iter, na.rm=TRUE)+1, xlab='Number of Dilations', ylab='#', verbose=FALSE)
}else if(hist=='sky'){
try({
if(!is.null(x$objects_redo)){
tempsky = image[x$objects_redo==0]
}else{
tempsky = image[x$objects==0]
}
tempsky = tempsky[tempsky > -8 & tempsky < 8 & !is.na(tempsky)]
if(!is.null(x$skyLL)){stat_LL = signif(x$skyLL,3)}else{stat_LL = NA}
if(!is.null(x$skyChiSq)){stat_ChiSq = signif(x$skyChiSq,3)}else{stat_ChiSq = NA}
magplot(function(x){dnorm(x,mean=0, sd=1)}, grid=TRUE, xlim=c(-6,6), xlab='(image - sky) / skyRMS', ylab='PDF', log='y', ylim=c(1e-8,0.5), lty=2, type='l', col='green4')
lines(density(tempsky[is.finite(tempsky)], bw=0.1, na.rm=TRUE, from=-8, to=8), col='black')
if(stat_sd_sky > 0){
lines(density((x$sky[x$objects==1] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='red')
lines(density((x$sky[x$objects==0] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='blue')
legend('bottom', legend=c('Sky pixels','Normal Dist','Sky in objects','Sky in sky'), lty=c(1,2,1,1), col=c('black','green4', 'red','blue'))
legend('topright',legend=c(paste0('Cor RMS: ',stat_cor_sky_skyRMS), paste0('Cor Image: ',stat_cor_sky_image)), bty='n', text.col='black')
}else{
legend('bottom', legend=c('Sky pixels','Normal Dist'), lty=c(1,2), col=c('black','green4'))
}
legend('topleft',legend=c(paste0('LL: ',stat_LL), paste0('Chi-Sq: ',stat_ChiSq)), bty='n', text.col='black')
})
}else{stop('Not a recognised hist type! Must be iters / sky.')}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 8 size v mag
if(logR50){
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(min(x$segstats$R50, 0.3, na.rm = TRUE), max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE, log='y')
}else{
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(0, max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE)
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 9 SB v rel_err
fluxrat=x$segstats$flux/x$segstats$flux_err
magplot(x$segstats$SB_N90, fluxrat, pch='.', col=hsv(alpha=0.5), ylim=c(0.5,max(fluxrat, 1, na.rm=TRUE)), cex=3, xlab='SB90 / mag/asec-sq', ylab='Flux/Flux-Error', grid=TRUE, log='y')
}else{
par(mar=c(3.5,3.5,0.5,0.5)) #plot 1 image
magimage(image, stretchscale=stretchscale, locut=-maximg, hicut=maximg, range=c(-1,1), type='num', zlim=c(-1,1), col=cmap)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 2 seg
# magimage(x$segim, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE)
# if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
profoundSegimPlot(profound=x)
#abline(v=c(0,dim(image)[1]))
#abline(h=c(0,dim(image)[2]))
par(mar=c(3.5,3.5,0.5,0.5)) #plot 3 dilate
magimage(image)
magimage(segdiff, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE, add=TRUE)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 4 mag hist
temphist=maghist(x$segstats$mag, log='y', scale=(2*dmag), breaks=seq(floor(min(x$segstats$mag, na.rm = TRUE)), ceiling(max(x$segstats$mag, na.rm = TRUE)),by=0.5), xlab='mag', ylab=paste('#/d',dmag,'mag',sep=''), grid=TRUE, verbose=FALSE)
ymax=log10(max(temphist$counts,na.rm = T))
xmax=temphist$mids[which.max(temphist$counts)]
abline(ymax - xmax*0.4, 0.4, col='red')
abline(v=xmax+0.25, col='red')
axis(side=1, at=xmax+0.25, labels=xmax+0.25, tick=FALSE, line=-1, col.axis='red')
legend('topleft', legend=paste('N:',length(x$segstats$mag)), bty='n', text.col='black')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 5 sky
magimage(x$sky-median(x$sky,na.rm=TRUE), qdiff=TRUE)
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky',paste0('Mean: ',stat_mean_sky),paste0('SD: ',stat_sd_sky)), bty='n', text.col='magenta')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 6 skyRMS
magimage(x$skyRMS)
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'orange'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky RMS',paste0('Mean: ',stat_mean_skyRMS),paste0('SD: ',stat_sd_skyRMS)), bty='n', text.col='magenta')
#plot 7 sky den
if(hist=='iters'){
maghist(x$segstats$iter, breaks=seq(-0.5,max(x$segstats$iter, na.rm=TRUE)+0.5,by=1), majorn=max(x$segstats$iter, na.rm=TRUE)+1, xlab='Number of Dilations', ylab='#', verbose=FALSE)
}else if(hist=='sky'){
try({
if(!is.null(x$objects_redo)){
tempsky = image[x$objects_redo==0]
}else{
tempsky = image[x$objects==0]
}
tempsky = tempsky[tempsky > -8 & tempsky < 8 & !is.na(tempsky)]
if(!is.null(x$skyLL)){stat_LL = signif(x$skyLL,3)}else{stat_LL = NA}
if(!is.null(x$skyChiSq)){stat_ChiSq = signif(x$skyChiSq,3)}else{stat_ChiSq = NA}
magplot(function(x){dnorm(x,mean=0, sd=1)}, grid=TRUE, xlim=c(-6,6), xlab='(image - sky) / skyRMS', ylab='PDF', log='y', ylim=c(1e-8,0.5), lty=2, type='l', col='green4')
lines(density(tempsky[is.finite(tempsky)], bw=0.1, na.rm=TRUE, from=-8, to=8), col='black')
if(stat_sd_sky > 0){
lines(density((x$sky[x$objects==1] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='red')
lines(density((x$sky[x$objects==0] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='blue')
legend('bottom', legend=c('Sky pixels','Normal Dist','Sky in objects','Sky in sky'), lty=c(1,2,1,1), col=c('black','green4', 'red','blue'))
legend('topright',legend=c(paste0('Cor RMS: ',stat_cor_sky_skyRMS), paste0('Cor Image: ',stat_cor_sky_image)), bty='n', text.col='black')
}else{
legend('bottom', legend=c('Sky pixels','Normal Dist'), lty=c(1,2), col=c('black','green4'))
}
legend('topleft',legend=c(paste0('LL: ',stat_LL), paste0('Chi-Sq: ',stat_ChiSq)), bty='n', text.col='black')
})
}else{stop('Not a recognised hist type! Must be iters / sky.')}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 8 size v mag
if(logR50){
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(min(x$segstats$R50, 0.3, na.rm = TRUE), max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE, log='y')
}else{
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(0, max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE)
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 9 SB v rel_err
fluxrat=x$segstats$flux/x$segstats$flux_err
magplot(x$segstats$SB_N90, fluxrat, pch='.', col=hsv(alpha=0.5), ylim=c(0.5,max(fluxrat, 1, na.rm=TRUE)), cex=3, xlab='SB90 / mag/pix-sq', ylab='Flux/Flux-Error', grid=TRUE, log='y')
}
})
}
.selectCoG=function(diffmat, threshold=1.05){
IDmat=matrix(rep(1:dim(diffmat)[2],each=dim(diffmat)[1]),nrow=dim(diffmat)[1])
logmat=diffmat>1 & diffmat<threshold
IDfin=IDmat
IDfin[logmat==FALSE]=NA
NegFlux=which(diffmat<threshold^0.2,arr.ind=TRUE)
if(length(NegFlux)>0){
NegFlux[,2]=NegFlux[,2]-1
IDfin[NegFlux]=IDmat[NegFlux]
IDfin[NegFlux[NegFlux[,2]==0,1],1]=0
}
tempout=suppressWarnings(apply(IDfin,1,min,na.rm=TRUE))
tempout[is.infinite(tempout)]=dim(diffmat)[2]
tempout+1
}
asgr/ProFound documentation built on Feb. 10, 2024, 9:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .selectCoG plot.profound profoundProFound
Documented in plot.profound profoundProFound
profoundProFound=function(image=NULL, segim=NULL, objects=NULL, mask=NULL, skycut=1.5, pixcut=3,
tolerance=4, ext=2, reltol=0, cliptol=Inf, sigma=1, smooth=TRUE,
SBlim=NULL, SBdilate=2, SBN100=100, size=5, shape='disc', iters=6,
threshold=1.05, magzero=0, gain=NULL, pixscale=1, sky=NULL, skyRMS=NULL,
redosegim=FALSE, redosky=TRUE, redoskysize=21, box=c(100,100), grid=box,
skygrid_type = 'new', type='bicubic', skytype='median', skyRMStype='quanlo', roughpedestal=FALSE,
sigmasel=1, skypixmin=prod(box)/2, boxadd=box/2, boxiters=0, conviters=100, iterskyloc=FALSE,
deblend=FALSE, df=3, radtrunc=2, iterative=FALSE, doChiSq=FALSE, doclip=TRUE,
shiftloc = FALSE, paddim = TRUE, keyvalues=NULL, verbose=FALSE,
plot=FALSE, stats=TRUE, rotstats=FALSE, boundstats=FALSE,
nearstats=boundstats, groupstats=boundstats, group=NULL,
groupby='segim_orig', offset=1, haralickstats=FALSE, sortcol="segID",
decreasing=FALSE, lowmemory=FALSE, keepim=TRUE, watershed='ProFound',
pixelcov=FALSE, deblendtype='fit', psf=NULL, fluxweight='sum',
convtype = 'brute', convmode = 'extended', fluxtype='Raw',
app_diam=1, Ndeblendlim=Inf, static_photom = FALSE, rem_mask=FALSE, ...){
if(verbose){message('Running ProFound:')}
timestart=proc.time()[3]
call=match.call(expand.dots=TRUE)
if(length(box)==1){
box=rep(box,2)
if(missing(grid)){grid=box}
if(missing(boxadd)){boxadd=box/2}
if(missing(skypixmin)){skypixmin=prod(box)/2}
}
if(length(grid)==1){
grid=rep(grid,2)
}
if(length(boxadd)==1){
boxadd=rep(boxadd,2)
}
fluxtype = tolower(fluxtype)
if(fluxtype=='raw' | fluxtype=='adu' | fluxtype=='adus'){
if(verbose){message('Using raw flux units')}
fluxscale=1
}else if (fluxtype=='jansky'){
if(verbose){message('Using Jansky flux units (WARNING: magzero must take system to AB)')}
fluxscale=10^(-0.4*(magzero-8.9))
}else if (fluxtype=='microjansky'){
if(verbose){message('Using Micro-Jansky flux units (WARNING: magzero must take system to AB)')}
fluxscale=10^(-0.4*(magzero-23.9))
}else{
stop('fluxtype must be Jansky / Microjansky / Raw!')
}
#Split out image and keyvalues parts of input:
if(!is.null(image)){
if(inherits(image, 'Rfits_image')){
keyvalues = image$keyvalues
image = image$imDat
}else if(inherits(image, 'Rfits_pointer')){
keyvalues = image$keyvalues
image = image[,]$imDat
}else if(inherits(image, 'matrix')){
'Do nothing'
}else{
stop('As of ProFound v1.21.0 only Rfits_image FITS inputs are allowed. Please install from GitHub asgr/Rfits')
}
}else{
stop('Missing image - this is a required input!')
}
if(!is.null(keyvalues)){
if(!inherits(keyvalues, 'Rfits_keylist')){
if(is.list(keyvalues)){
class(keyvalues) = 'Rfits_keylist'
}else{
stop('keyvalues is the wrong format- should be a list!')
}
}
}
if(box[1] > ceiling(dim(image)[1]/3)){
box[1] = ceiling(dim(image)[1]/3)
message('dim(image)[1]/box[1] must be >=3, box[1] modified to ',box[1])
}
if(box[2] > ceiling(dim(image)[2]/3)){
box[2] = ceiling(dim(image)[2]/3)
message('dim(image)[2]/box[2] must be >=3, box[2] modified to ',box[2])
}
if(grid[1] > ceiling(dim(image)[1]/3)){
grid[1] = ceiling(dim(image)[1]/3)
message('dim(image)[1]/grid[1] must be >=3, grid[1] modified to ',grid[1])
}
if(grid[2] > ceiling(dim(image)[2]/3)){
grid[2] = ceiling(dim(image)[2]/3)
message('dim(image)[2]/grid[2] must be >=3, grid[2] modified to ',grid[2])
}
if(verbose){message(paste('Supplied image is',dim(image)[1],'x',dim(image)[2],'pixels'))}
#Treat image NAs as masked regions:
badpix=NULL
if(!is.null(mask)){
mask=mask*1L #Looks silly, but this ensures a logical mask becomes integer.
if(length(mask)==1 & !is.na(mask[1])){
maskflag=mask
mask=matrix(0L,dim(image)[1],dim(image)[2])
mask[image==maskflag]=1L
}
if(anyNA(image)){
badpix=which(is.na(image))
mask[badpix]=1L
image[badpix]=0
}
if(is.numeric(mask)){
mask=as.integer(mask)
attributes(mask)$dim = dim(image)
}
}else{
if(anyNA(image)){
mask=matrix(0L,dim(image)[1],dim(image)[2])
badpix=which(is.na(image))
mask[badpix]=1L
image[badpix]=0
}
}
#if(!is.null(segim) & !is.null(mask)){
# segim=segim*(1-mask) #I don't think we actually need this
#}
#Get the pixel scale, if possible and not provided:
if(missing(pixscale) & !is.null(keyvalues)){
pixscale = pixscale(keyvalues)
if(verbose){message(paste('Extracted pixel scale from keyvalues provided:',signif(pixscale,4),'asec/pixel'))}
}else{
if(verbose){message(paste('Using suggested pixel scale:',signif(pixscale,4),'asec/pixel'))}
}
imarea = prod(dim(image))*pixscale^2/(3600^2)
if(verbose){message(paste('Supplied image is',signif(dim(image)[1]*pixscale/60,4),'x',signif(dim(image)[2]*pixscale/60,4),'amin, ', signif(imarea,4),'deg-sq'))}
if(is.null(objects)){
if(!is.null(segim)){
objects = matrix(0L,dim(segim)[1],dim(segim)[2])
objects[] = as.logical(segim)
}
}#else{
# objects = objects*1L #Looks silly, but this ensures a logical mask becomes integer.
# }
if(static_photom & !is.null(segim)){
#set all the flags to make static photom
sky = 0
iters = 0
redosky = FALSE
redosegim = FALSE
}
#Check for user provided sky, and compute if missing:
hassky = !is.null(sky)
hasskyRMS = !is.null(skyRMS)
if((hassky==FALSE | hasskyRMS==FALSE) & is.null(segim)){
if(hassky==FALSE){
sky = 0
}
if(verbose){message(paste('Making initial sky map -',round(proc.time()[3]-timestart,3),'sec'))}
roughsky = profoundMakeSkyGrid(image=image, objects=objects, mask=mask, sky=sky, box=box,
grid=grid, skygrid_type=skygrid_type, type=type,
skytype=skytype, skyRMStype=skyRMStype, sigmasel=sigmasel,
skypixmin=skypixmin, boxadd=boxadd, boxiters=0, conviters=conviters,
doclip=doclip, shiftloc=shiftloc, paddim=paddim, rem_mask=rem_mask)
if(roughpedestal){
roughsky$sky = median(roughsky$sky,na.rm=TRUE)
roughsky$skyRMS = median(roughsky$skyRMS,na.rm=TRUE)
}
if(hassky==FALSE){
sky = roughsky$sky
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
}
if(hasskyRMS==FALSE){
skyRMS = roughsky$skyRMS
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}
rm(roughsky)
}else{
if(verbose){message("Skipping making initial sky map - User provided sky and sky RMS, or user provided segim")}
}
#Make the initial segmentation map, if not provided.
if(is.null(segim)){
if(verbose){message(paste('Making initial segmentation image -',round(proc.time()[3]-timestart,3),'sec'))}
segim = profoundMakeSegim(image=image, objects=objects, mask=mask, tolerance=tolerance,
ext=ext, reltol=reltol, cliptol=cliptol, sigma=sigma,
smooth=smooth, pixcut=pixcut, skycut=skycut, SBlim=SBlim,
sky=sky, skyRMS=skyRMS, magzero=magzero, pixscale=pixscale,
verbose=verbose, watershed=watershed, plot=FALSE, stats=FALSE)
objects = segim$objects
segim = segim$segim
}else{
redosegim = FALSE
if(verbose){message("Skipping making an initial segmentation image - User provided segim")}
}
if(any(segim>0)){
if((hassky==FALSE | hasskyRMS==FALSE)){
if(redosky){
if(verbose){message(paste('Doing initial aggressive dilation -',round(proc.time()[3]-timestart,3),'sec'))}
objects_redo=profoundMakeSegimDilate(segim=objects, size=redoskysize, shape=shape, sky=sky, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$objects
}else{
objects_redo=objects
}
if(verbose){message(paste('Making better sky map -',round(proc.time()[3]-timestart,3),'sec'))}
if(hasskyRMS==FALSE){rm(skyRMS)}
bettersky = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky,
box=box, skygrid_type=skygrid_type, grid=grid,
type=type, skytype=skytype, skyRMStype=skyRMStype,
sigmasel=sigmasel, skypixmin=skypixmin, boxadd=boxadd,
boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)
if(hassky==FALSE){
sky = bettersky$sky
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
}
if(hasskyRMS==FALSE){
skyRMS = bettersky$skyRMS
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}
rm(bettersky)
if(redosegim){
if(verbose){message(paste('Making better segmentation image -',round(proc.time()[3]-timestart,3),'sec'))}
imagescale = (image-sky)/skyRMS
imagescale[!is.finite(imagescale)] = 0
if(!is.null(SBlim) & !missing(magzero)){
imagescale[imagescale<skycut | sky<profoundSB2Flux(SBlim, magzero, pixscale)]=0
}else{
imagescale[imagescale<skycut] = 0
}
if(!is.null(mask)){
imagescale[mask!=0] = 0
}
segim[imagescale==0] = 0
objects[segim==0] = 0
rm(imagescale)
}
}else{
if(verbose){message("Skipping making better sky map - User provided sky and sky RMS")}
}
if(iters>0 | iterskyloc){
if(verbose){message(paste('Calculating initial segstats -',round(proc.time()[3]-timestart,3),'sec'))}
# if(length(sky)>1){
# skystats = .profoundFluxCalcMin(image = sky - median(sky,na.rm=TRUE), segim=segim, mask=mask) #run on sky
# skystats = skystats$flux / skystats$N100 #get per pixel mean flux per segment
# #skymed = median(skystats, na.rm=TRUE) #median per pixel mean flux per segment for the whole image
# }else{
# skystats = 0 #specified
# #skymed = sky #specified
# }
image_sky = image - sky
segstats = .profoundFluxCalcMin(image=image_sky, segim=segim, mask=mask) #run on initial image
#segstats$flux = segstats$flux - (skystats * segstats$N100) #remove the local sky component
origfrac = segstats$flux #set to initial fluxes
segim_orig = segim
expand_segID = segstats[,'segID']
SBlast = rep(Inf,length(expand_segID))
selseg = rep(0,length(expand_segID))
if(is.null(SBdilate)){
SBdilate = Inf
}else{
skyRMSstats = .profoundFluxCalcMin(image=skyRMS, segim=segim, mask=mask) #run on sky
skyRMSstats = skyRMSstats$flux/skyRMSstats$N100 #get per pixel mean flux per segment
SBdilate = skyRMSstats * 10^(-0.4*SBdilate)
SBdilate[!is.finite(SBdilate)]=Inf
}
if(iters>0){
if(verbose){message('Doing dilations:')}
for(i in 1:(iters)){
if(verbose){message(paste('Iteration',i,'of',iters,'-',round(proc.time()[3]-timestart,3),'sec'))}
segim_new = profoundDilate(segim=segim, size=size, shape=shape, expand=expand_segID, iters=1)
#segim_new = profoundMakeSegimDilate(segim=segim, expand=expand_segID, size=size, shape=shape, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$segim #dilate
segstats_new = .profoundFluxCalcMin(image=image_sky, segim=segim_new, mask=mask) #run on image with dilated segments
N100diff = (segstats_new$N100 - segstats$N100)
SBnew = (segstats_new$flux - segstats$flux) / N100diff #calculate the surface brightness of the new grown anulus only
fluxgrowthratio = segstats_new$flux / segstats$flux #calculate flux growth ratio
#skyfrac = abs((skystats * N100diff) / (segstats_new$flux - segstats$flux)) #estimate how much of the flux growth might be coming from unusual local sky
expand_segID = segstats[which((fluxgrowthratio > threshold | (SBnew > SBdilate & N100diff>SBN100)) & segstats_new$flux>0 & SBnew < SBlast & selseg==(i-1)),'segID']
expand_segID = expand_segID[is.finite(expand_segID)] #safety first
if(length(expand_segID)==0){break}
updateID = which(segstats$segID %in% expand_segID)
selseg[updateID] = i #iteration number for segments flagged for growth
segstats[updateID,] = segstats_new[updateID,] #update only the segstats for things that are growing, the rest can be ignored
SBlast = SBnew
if('fastmatch' %in% .packages()){ #dilate segments that pass tests
selpix = which(fastmatch::fmatch(segim_new, expand_segID, nomatch = 0L) > 0)
}else{
selpix = which(segim_new %in% expand_segID)
}
segim[selpix] = segim_new[selpix]
}
}
origfrac = origfrac / segstats$flux
if(iterskyloc){
segim_skyloc = profoundMakeSegimDilate(segim=segim, size=size, shape=shape, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$segim
segstats_sky = .profoundFluxCalcMin(image=image, segim=segim_skyloc-segim, mask=mask)
#segstats$flux = segstats$flux + (skystats * segstats$N100) #add back the local sky component
skyseg_temp = segstats_sky$flux / segstats_sky$N100
skyseg_temp[!is.finite(skyseg_temp)] = 0
skyseg_mean = rep(NA,length(segstats$segID))
skyseg_mean[match(segstats_sky$segID,segstats$segID)] = skyseg_temp
}else{
skyseg_mean = NA
}
objects = matrix(0L,dim(segim)[1],dim(segim)[2])
objects[] = as.logical(segim)
}else{
if(verbose){message('Iters set to 0 - keeping segim un-dilated')}
segim_orig = segim
selseg = 0
origfrac = 1
skyseg_mean = NA
}
objects_redo = profoundMakeSegimDilate(segim=objects, mask=mask, size=redoskysize, shape=shape, sky=sky, verbose=verbose, plot=FALSE, stats=FALSE, rotstats=FALSE)$objects
if(redosky){
if(redoskysize %% 2 == 0){redoskysize=redoskysize+1}
if(verbose){message(paste('Doing final aggressive dilation -',round(proc.time()[3]-timestart,3),'sec'))}
if(verbose){message(paste('Making final sky map -',round(proc.time()[3]-timestart,3),'sec'))}
rm(skyRMS)
sky = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type=skygrid_type, grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)$sky
if(verbose){message(paste('Making final sky RMS map -',round(proc.time()[3]-timestart,3),'sec'))}
skyRMS = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type=skygrid_type, grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doclip=doclip, shiftloc=shiftloc,
paddim=paddim, rem_mask=rem_mask)$skyRMS
if(verbose){message(' - Sky statistics :')}
if(verbose){print(summary(as.numeric(sky[!is.na(sky)])))}
if(verbose){message(' - Sky-RMS statistics :')}
if(verbose){print(summary(as.numeric(skyRMS[!is.na(skyRMS)])))}
}else{
if(verbose){message("Skipping making final sky map - redosky set to FALSE")}
}
if(doChiSq){
if(verbose){message(paste('Making sky LL ChiSq map -',round(proc.time()[3]-timestart,3),'sec'))}
skyChiSqMap = profoundMakeSkyGrid(image=image, objects=objects_redo, mask=mask, sky=sky, box=box,
skygrid_type='old', grid=grid, type=type, skytype=skytype,
skyRMStype=skyRMStype, sigmasel=sigmasel, skypixmin=skypixmin,
boxadd=boxadd, boxiters=boxiters, conviters=conviters, doChiSq=TRUE,
doclip=doclip, shiftloc=shiftloc, paddim=paddim, rem_mask=rem_mask)$skyChiSq
}else{
skyChiSqMap = NA
}
Norig = tabulate(segim_orig)
if(is.function(pixelcov)){
cor_err_func = pixelcov
}else{
if(pixelcov){
if(verbose){message(paste('Calculating pixel covariance -',round(proc.time()[3]-timestart,3),'sec'))}
if(is.null(profoundMakeSkyGrid)){
cor_err_func = profoundPixelCorrelation(image=image, objects=objects, mask=mask, sky=sky, skyRMS=skyRMS,
fft=FALSE, lag=apply(expand.grid(c(1,2,4),c(1,10,100,1000,1e4)),MARGIN=1,FUN=prod))$cor_err_func
}else{
cor_err_func = profoundPixelCorrelation(image=image, objects=objects_redo, mask=mask, sky=sky, skyRMS=skyRMS,
fft=FALSE, lag=apply(expand.grid(c(1,2,4),c(1,10,100,1000,1e4)),MARGIN=1,FUN=prod))$cor_err_func
}
}else{
cor_err_func = NULL
}
}
if(lowmemory){
image = image - sky
sky = 0
skyRMS = 0
segim_orig = NULL
objects = NULL
objects_redo = NULL
}
if(stats & !is.null(image)){
if(verbose){message(paste('Calculating final segstats for',length(which(tabulate(segim)>0)),'objects -',round(proc.time()[3]-timestart,3),'sec'))}
if(verbose){message(paste(' - magzero =', signif(magzero,4)))}
if(verbose){
if(is.null(gain)){
message(paste(' - gain = NULL (ignored)'))
}else{
message(paste(' - gain =', signif(gain,4)))
}
}
if(verbose){message(paste(' - pixscale =', signif(pixscale,4)))}
if(verbose){message(paste(' - rotstats =', rotstats))}
if(verbose){message(paste(' - boundstats =', boundstats))}
segstats = profoundSegimStats(image=image, segim=segim, mask=mask, sky=sky, skyRMS=skyRMS,
magzero=magzero, gain=gain, pixscale=pixscale, keyvalues=keyvalues,
sortcol=sortcol, decreasing=decreasing, rotstats=rotstats,
boundstats=boundstats, offset=offset, cor_err_func=cor_err_func,
app_diam=app_diam)
segstats = cbind(segstats, iter=selseg, origfrac=origfrac, Norig=Norig[segstats$segID], skyseg_mean=skyseg_mean)
segstats = cbind(segstats, flag_keep=segstats$origfrac >= median(segstats$origfrac[segstats$iter==iters]) | segstats$iter<iters)
}else{
if(verbose){message("Skipping segmentation statistics - segstats set to FALSE")}
segstats = NULL
}
if(nearstats){
near = profoundSegimNear(segim=segim, offset=offset)
}else{
near = NULL
}
if(deblend){
groupstats = TRUE
}
if(groupstats){
if(verbose){message(paste(' - groupstats = TRUE - ',round(proc.time()[3]-timestart,3),'sec'))}
if(groupby=='segim'){
if(is.null(group)){
group=profoundSegimGroup(segim)
}
}else if(groupby=='segim_orig'){
if(is.null(group)){
#message(round(proc.time()[3]-timestart,3))
group = profoundSegimGroup(segim_orig)
if(any(group$groupsegID$Ngroup>1)){
#message(round(proc.time()[3]-timestart,3))
group$groupim = profoundSegimKeep(segim=segim, segID_merge=group$groupsegID[group$groupsegID$Ngroup>1,'segID'])
#message(round(proc.time()[3]-timestart,3))
group$groupsegID$Npix = tabulate(group$groupim)[group$groupsegID$groupID]
}
}
}else{
stop('Non legal groupby option, must be segim or segim_orig!')
}
if(stats & !is.null(image) & !is.null(group)){
groupstats = profoundSegimStats(image=image, segim=group$groupim, mask=mask,
sky=sky, skyRMS=skyRMS, magzero=magzero, gain=gain,
pixscale=pixscale, keyvalues=keyvalues, sortcol=sortcol,
decreasing=decreasing, rotstats=rotstats, boundstats=boundstats,
offset=offset, cor_err_func=cor_err_func, app_diam=app_diam)
colnames(groupstats)[1] = 'groupID'
}else{
groupstats = NULL
}
}else{
if(verbose){message(' - groupstats = FALSE')}
group = NULL
groupstats = NULL
}
if(deblend & stats & !is.null(image) & any(group$groupsegID$Ngroup>1)){
if(verbose){message(paste(' - deblend = TRUE - ',round(proc.time()[3]-timestart,3),'sec'))}
tempblend = profoundFluxDeblend(image=image-sky, segim=segim, segstats=segstats,
groupim=group$groupim, groupsegID=group$groupsegID,
magzero=magzero, df=df, radtrunc=radtrunc, iterative=iterative,
doallstats=TRUE, deblendtype=deblendtype, psf=psf,
fluxweight=fluxweight, convtype=convtype, convmode=convmode,
Ndeblendlim = Ndeblendlim)
if(!is.null(tempblend)){
segstats = cbind(segstats,tempblend[,-2])
}
}else{
if(verbose){message(' - deblend = FALSE')}
}
if(haralickstats){
if(requireNamespace("EBImage", quietly = TRUE)){
scale = 10^(0.4*(30-magzero))
temphara = (image-sky)*scale
if(!is.null(mask)){
temphara[mask!=0] = 0
}
temphara[!is.finite(temphara)] = 0
haralick = as.data.frame(EBImage::computeFeatures.haralick(segim,temphara))
haralick = haralick[segstats$segID,]
}else{
if(verbose){
message('The EBImage package is needed to compute Haralick statistics.')
haralick = NULL
}
}
}else{
haralick = NULL
}
if(plot){
if(verbose){message(paste('Plotting segments -',round(proc.time()[3]-timestart,3),'sec'))}
if(any(is.finite(sky))){
profoundSegimPlot(image=image-sky, segim=segim, mask=mask, keyvalues=keyvalues, ...)
}else{
profoundSegimPlot(image=image, segim=segim, mask=mask, keyvalues=keyvalues, ...)
}
}else{
if(verbose){message("Skipping segmentation plot - plot = FALSE")}
}
if(is.null(SBlim)){
SBlim = NULL
}else if(is.numeric(SBlim)){
SBlimtemp = profoundFlux2SB(flux=skyRMS*skycut, magzero=magzero, pixscale=pixscale)
SBlimtemp = matrix(SBlimtemp,dim(skyRMS)[1],dim(skyRMS)[2])
SBlimtemp[which(SBlimtemp>SBlim)] = SBlim
SBlim = SBlimtemp
}else if(SBlim[1]=='get' & skycut> -Inf){
SBlim = profoundFlux2SB(flux=skyRMS*skycut, magzero=magzero, pixscale=pixscale)
}
if(is.null(keyvalues)){keyvalues = NULL}
if(keepim==FALSE){image = NULL; mask = NULL}
if(is.null(mask)){mask = NULL}
if(!is.null(badpix)){image[badpix] = NA}
row.names(segstats) = NULL
segstats[,grep('flux',colnames(segstats))] = fluxscale*segstats[,grep('flux',colnames(segstats))]
if(!is.null(groupstats)){
groupstats[,grep('flux',colnames(groupstats))] = fluxscale*groupstats[,grep('flux',colnames(groupstats))]
}
if(stats){
cutsky = (image - sky) / skyRMS
cutsky[!is.finite(cutsky)] = NA
if(!is.null(mask)){
cutsky[mask>0] = NA
}
if(!is.null(objects_redo)){
cutsky[objects_redo == 1] = NA
}else{
cutsky[objects == 1] = NA
}
cutsky = cutsky[which(cutsky<=0)] #Only consider sky pixels below the mode
if(length(cutsky) > 0){
df = length(cutsky) - 1
skyChiSq = sum(cutsky^2, na.rm =TRUE)
skyLL = dchisq(skyChiSq, df=df, log=TRUE)
skyChiSq = skyChiSq / df
}else{
skyLL = NA
skyChiSq = NA
}
}else{
skyLL = NULL
skyChiSq = NULL
}
if(verbose){message(paste('ProFound is finished! -',round(proc.time()[3]-timestart,3),'sec'))}
output = list(segim=segim, segim_orig=segim_orig, objects=objects, objects_redo=objects_redo,
sky=sky, skyRMS=skyRMS, image=image, mask=mask, segstats=segstats,
Nseg=dim(segstats)[1], near=near, group=group, groupstats=groupstats,
haralick=haralick, keyvalues=keyvalues, SBlim=SBlim, magzero=magzero, dim=dim(segim),
pixscale=pixscale, imarea=imarea, skyLL=skyLL, skyChiSq=skyChiSq, skyChiSqMap=skyChiSqMap,
gain=gain, call=call, date=date(), time=proc.time()[3]-timestart, ProFound.version=packageVersion('ProFound'),
R.version=R.version)
}else{
if(is.null(keyvalues)){keyvalues = NULL}
if(keepim==FALSE){image = NULL; mask = NULL}
if(is.null(mask)){mask = NULL}
if(!is.null(badpix)){image[badpix] = NA}
if(verbose){message('No objects in segmentation map - skipping dilations and CoG')}
if(verbose){message(paste('ProFound is finished! -',round(proc.time()[3]-timestart,3),'sec'))}
output = list(segim=NULL, segim_orig=NULL, objects=NULL, objects_redo=NULL, sky=sky,
skyRMS=skyRMS, image=image, mask=mask, segstats=NULL, Nseg=0, near=NULL,
group=NULL, groupstats=NULL, haralick=NULL, keyvalues=keyvalues, SBlim=NULL,
magzero=magzero, dim=dim(segim), pixscale=pixscale, imarea=imarea, skyLL=NULL, skyChiSq=NULL,
skyChiSqMap=NULL, gain=gain, call=call, date=date(), time=proc.time()[3]-timestart,
ProFound.version=packageVersion('ProFound'), R.version=R.version)
}
class(output) = 'profound'
return(invisible(output))
}
plot.profound=function(x, logR50=TRUE, dmag=0.5, hist='sky', ...){
suppressWarnings({
if(!inherits(x, 'profound')){
stop('Object class is not of type profound!')
}
if(is.null(x$image)){
stop('Missing image!')
}
if(is.null(x$segim)){
stop('Missing segmentation map!')
}
if(is.null(x$sky)){
x$sky=matrix(0, x$dim[1], x$dim[2])
}
if(length(x$sky)==1){
x$sky=matrix(x$sky, x$dim[1], x$dim[2])
}
if(is.null(x$skyRMS)){
x$skyRMS=matrix(1, x$dim[1], x$dim[2])
}
if(length(x$skyRMS)==1){
x$skyRMS=matrix(x$skyRMS, x$dim[1], x$dim[2])
}
segdiff = x$segim - x$segim_orig
segdiff[segdiff<0]=0
if(!is.null(x$mask)){
x$image[x$mask>0] = NA
x$mask[is.na(x$image)] = 1
}
image = (x$image-x$sky)/x$skyRMS
image[x$skyRMS < 0] = NA
if(is.null(x$objects)){
x$objects = matrix(0L,dim(x$segim)[1],dim(x$segim)[2])
x$objects[] = as.logical(x$segim)
}
#Masked and unmasked stats:
if(!is.null(x$mask)){
stat_mean_sky = signif(mean(x$sky[x$mask==0], na.rm=TRUE),4)
stat_sd_sky = signif(sd(x$sky[x$mask==0], na.rm=TRUE),4)
stat_mean_skyRMS = signif(mean(x$skyRMS[x$mask==0], na.rm=TRUE),4)
stat_sd_skyRMS = signif(sd(x$skyRMS[x$mask==0], na.rm=TRUE),4)
stat_cor_sky_skyRMS = signif(cor(as.numeric(x$sky[x$mask==0]), as.numeric(x$skyRMS[x$mask==0])^2, use="pairwise.complete.obs"),4)
if(!is.null(x$objects_redo)){
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$mask==0 & x$objects_redo>0]), as.numeric(x$image[x$mask==0 & x$objects_redo>0] - as.numeric(x$sky[x$mask==0 & x$objects_redo>0])), use="pairwise.complete.obs"),4)
}else{
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$mask==0 & x$objects>0]), as.numeric(x$image[x$mask==0 & x$objects>0] - as.numeric(x$sky[x$mask==0 & x$objects>0])), use="pairwise.complete.obs"),4)
}
}else{
stat_mean_sky = signif(mean(x$sky, na.rm=TRUE),4)
stat_sd_sky = signif(sd(x$sky, na.rm=TRUE),4)
stat_mean_skyRMS = signif(mean(x$skyRMS, na.rm=TRUE),4)
stat_sd_skyRMS = signif(sd(x$skyRMS, na.rm=TRUE),4)
stat_cor_sky_skyRMS = signif(cor(as.numeric(x$sky), as.numeric(x$skyRMS)^2, use="pairwise.complete.obs"),4)
if(!is.null(x$objects_redo)){
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$objects_redo>0]), as.numeric(x$image[x$objects_redo>0] - as.numeric(x$sky[x$objects_redo>0])), use="pairwise.complete.obs"),4)
}else{
stat_cor_sky_image = signif(cor(as.numeric(x$sky[x$objects>0]), as.numeric(x$image[x$objects>0] - as.numeric(x$sky[x$objects>0])), use="pairwise.complete.obs"),4)
}
}
if(requireNamespace("imager", quietly = TRUE)){
outline = x$objects - as.matrix(imager::erode(imager::as.cimg(x$objects),imager::imfill(3,3,val=1)))
}else{
outline = NULL
}
cmap = hcl.colors(100, 'RdYlBu', rev=TRUE)
maximg = quantile(abs(image[is.finite(image)]), 0.995, na.rm=TRUE)
stretchscale = 1/median(abs(image), na.rm=TRUE)
layout(matrix(1:9, 3, byrow=TRUE))
if(!is.null(x$keyvalues)){
par(mar=c(3.5,3.5,0.5,0.5)) #plot 1 image
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=image, keyvalues=x$keyvalues, stretchscale=stretchscale, locut=-maximg, hicut=maximg, range=c(-1,1), type='num', zlim=c(-1,1), col=cmap)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 2 seg
profoundSegimPlot(profound=x)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 3 dilate
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=image, keyvalues=x$keyvalues)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
magimage(segdiff, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE, add=TRUE)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 4 mag hist
if(is.null(x$imarea)){
imarea=prod(x$dim)*x$pixscale^2/(3600^2)
}else{
imarea=x$imarea
}
temphist=maghist(x$segstats$mag, log='y', scale=(2*dmag)/x$imarea, breaks=seq(floor(min(x$segstats$mag, na.rm = TRUE)), ceiling(max(x$segstats$mag, na.rm = TRUE)),by=0.5),
xlab='mag', ylab=paste('#/deg-sq/d',dmag,'mag',sep=''), grid=TRUE, verbose=FALSE)
#magplot(temphist, log='y', xlab='mag', ylab=expression('#'/'deg-sq'/'dmag'), grid=TRUE)
ymax=log10(max(temphist$counts,na.rm = T))
xmax=temphist$mids[which.max(temphist$counts)]
abline(ymax - xmax*0.4, 0.4, col='red')
abline(v=xmax+0.25, col='red')
axis(side=1, at=xmax+0.25, labels=xmax+0.25, tick=FALSE, line=-1, col.axis='red')
legend('topleft', legend=paste('N:',length(x$segstats$mag)), bty='n', text.col='black')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 5 sky
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=x$sky-median(x$sky,na.rm=TRUE), keyvalues=x$keyvalues, qdiff=TRUE)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky',paste0('Mean: ',stat_mean_sky),paste0('SD: ',stat_sd_sky)), bty='n', text.col='magenta')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 6 skyRMS
if(requireNamespace("Rwcs", quietly = TRUE)){
Rwcs::Rwcs_image(image=x$skyRMS, keyvalues=x$keyvalues)
}else{
stop("The Rwcs package is need to process the keyvalues. Install from GitHub asgr/Rfits.")
}
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'orange'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky RMS',paste0('Mean: ',stat_mean_skyRMS),paste0('SD: ',stat_sd_skyRMS)), bty='n', text.col='magenta')
#plot 7 sky den
if(hist=='iters'){
maghist(x$segstats$iter, breaks=seq(-0.5,max(x$segstats$iter, na.rm=TRUE)+0.5,by=1), majorn=max(x$segstats$iter, na.rm=TRUE)+1, xlab='Number of Dilations', ylab='#', verbose=FALSE)
}else if(hist=='sky'){
try({
if(!is.null(x$objects_redo)){
tempsky = image[x$objects_redo==0]
}else{
tempsky = image[x$objects==0]
}
tempsky = tempsky[tempsky > -8 & tempsky < 8 & !is.na(tempsky)]
if(!is.null(x$skyLL)){stat_LL = signif(x$skyLL,3)}else{stat_LL = NA}
if(!is.null(x$skyChiSq)){stat_ChiSq = signif(x$skyChiSq,3)}else{stat_ChiSq = NA}
magplot(function(x){dnorm(x,mean=0, sd=1)}, grid=TRUE, xlim=c(-6,6), xlab='(image - sky) / skyRMS', ylab='PDF', log='y', ylim=c(1e-8,0.5), lty=2, type='l', col='green4')
lines(density(tempsky[is.finite(tempsky)], bw=0.1, na.rm=TRUE, from=-8, to=8), col='black')
if(stat_sd_sky > 0){
lines(density((x$sky[x$objects==1] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='red')
lines(density((x$sky[x$objects==0] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='blue')
legend('bottom', legend=c('Sky pixels','Normal Dist','Sky in objects','Sky in sky'), lty=c(1,2,1,1), col=c('black','green4', 'red','blue'))
legend('topright',legend=c(paste0('Cor RMS: ',stat_cor_sky_skyRMS), paste0('Cor Image: ',stat_cor_sky_image)), bty='n', text.col='black')
}else{
legend('bottom', legend=c('Sky pixels','Normal Dist'), lty=c(1,2), col=c('black','green4'))
}
legend('topleft',legend=c(paste0('LL: ',stat_LL), paste0('Chi-Sq: ',stat_ChiSq)), bty='n', text.col='black')
})
}else{stop('Not a recognised hist type! Must be iters / sky.')}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 8 size v mag
if(logR50){
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(min(x$segstats$R50, 0.3, na.rm = TRUE), max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE, log='y')
}else{
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(0, max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE)
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 9 SB v rel_err
fluxrat=x$segstats$flux/x$segstats$flux_err
magplot(x$segstats$SB_N90, fluxrat, pch='.', col=hsv(alpha=0.5), ylim=c(0.5,max(fluxrat, 1, na.rm=TRUE)), cex=3, xlab='SB90 / mag/asec-sq', ylab='Flux/Flux-Error', grid=TRUE, log='y')
}else{
par(mar=c(3.5,3.5,0.5,0.5)) #plot 1 image
magimage(image, stretchscale=stretchscale, locut=-maximg, hicut=maximg, range=c(-1,1), type='num', zlim=c(-1,1), col=cmap)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 2 seg
# magimage(x$segim, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE)
# if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
profoundSegimPlot(profound=x)
#abline(v=c(0,dim(image)[1]))
#abline(h=c(0,dim(image)[2]))
par(mar=c(3.5,3.5,0.5,0.5)) #plot 3 dilate
magimage(image)
magimage(segdiff, col=c(NA, rainbow(max(x$segim,na.rm=TRUE), end=2/3)), magmap=FALSE, add=TRUE)
if(!is.null(x$mask)){magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 4 mag hist
temphist=maghist(x$segstats$mag, log='y', scale=(2*dmag), breaks=seq(floor(min(x$segstats$mag, na.rm = TRUE)), ceiling(max(x$segstats$mag, na.rm = TRUE)),by=0.5), xlab='mag', ylab=paste('#/d',dmag,'mag',sep=''), grid=TRUE, verbose=FALSE)
ymax=log10(max(temphist$counts,na.rm = T))
xmax=temphist$mids[which.max(temphist$counts)]
abline(ymax - xmax*0.4, 0.4, col='red')
abline(v=xmax+0.25, col='red')
axis(side=1, at=xmax+0.25, labels=xmax+0.25, tick=FALSE, line=-1, col.axis='red')
legend('topleft', legend=paste('N:',length(x$segstats$mag)), bty='n', text.col='black')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 5 sky
magimage(x$sky-median(x$sky,na.rm=TRUE), qdiff=TRUE)
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'black'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky',paste0('Mean: ',stat_mean_sky),paste0('SD: ',stat_sd_sky)), bty='n', text.col='magenta')
par(mar=c(3.5,3.5,0.5,0.5)) #plot 6 skyRMS
magimage(x$skyRMS)
if(!is.null(x$mask)){
magimage(x$mask!=0, col=c(NA,hsv(alpha=0.2)), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
if(!is.null(outline)){
magimage(outline, col=c(NA,'orange'), add=TRUE, magmap=FALSE, zlim=c(0,1))
}
legend('topleft',legend=c('Sky RMS',paste0('Mean: ',stat_mean_skyRMS),paste0('SD: ',stat_sd_skyRMS)), bty='n', text.col='magenta')
#plot 7 sky den
if(hist=='iters'){
maghist(x$segstats$iter, breaks=seq(-0.5,max(x$segstats$iter, na.rm=TRUE)+0.5,by=1), majorn=max(x$segstats$iter, na.rm=TRUE)+1, xlab='Number of Dilations', ylab='#', verbose=FALSE)
}else if(hist=='sky'){
try({
if(!is.null(x$objects_redo)){
tempsky = image[x$objects_redo==0]
}else{
tempsky = image[x$objects==0]
}
tempsky = tempsky[tempsky > -8 & tempsky < 8 & !is.na(tempsky)]
if(!is.null(x$skyLL)){stat_LL = signif(x$skyLL,3)}else{stat_LL = NA}
if(!is.null(x$skyChiSq)){stat_ChiSq = signif(x$skyChiSq,3)}else{stat_ChiSq = NA}
magplot(function(x){dnorm(x,mean=0, sd=1)}, grid=TRUE, xlim=c(-6,6), xlab='(image - sky) / skyRMS', ylab='PDF', log='y', ylim=c(1e-8,0.5), lty=2, type='l', col='green4')
lines(density(tempsky[is.finite(tempsky)], bw=0.1, na.rm=TRUE, from=-8, to=8), col='black')
if(stat_sd_sky > 0){
lines(density((x$sky[x$objects==1] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='red')
lines(density((x$sky[x$objects==0] - stat_mean_sky) / stat_sd_sky, na.rm=TRUE, from=-8, to=8), col='blue')
legend('bottom', legend=c('Sky pixels','Normal Dist','Sky in objects','Sky in sky'), lty=c(1,2,1,1), col=c('black','green4', 'red','blue'))
legend('topright',legend=c(paste0('Cor RMS: ',stat_cor_sky_skyRMS), paste0('Cor Image: ',stat_cor_sky_image)), bty='n', text.col='black')
}else{
legend('bottom', legend=c('Sky pixels','Normal Dist'), lty=c(1,2), col=c('black','green4'))
}
legend('topleft',legend=c(paste0('LL: ',stat_LL), paste0('Chi-Sq: ',stat_ChiSq)), bty='n', text.col='black')
})
}else{stop('Not a recognised hist type! Must be iters / sky.')}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 8 size v mag
if(logR50){
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(min(x$segstats$R50, 0.3, na.rm = TRUE), max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE, log='y')
}else{
magplot(x$segstats$mag, x$segstats$R50, pch='.', col=hsv(alpha=0.5), ylim=c(0, max(x$segstats$R50, 1, na.rm = TRUE)), cex=3, xlab='mag', ylab='R50 / asec', grid=TRUE)
}
par(mar=c(3.5,3.5,0.5,0.5)) #plot 9 SB v rel_err
fluxrat=x$segstats$flux/x$segstats$flux_err
magplot(x$segstats$SB_N90, fluxrat, pch='.', col=hsv(alpha=0.5), ylim=c(0.5,max(fluxrat, 1, na.rm=TRUE)), cex=3, xlab='SB90 / mag/pix-sq', ylab='Flux/Flux-Error', grid=TRUE, log='y')
}
})
}
.selectCoG=function(diffmat, threshold=1.05){
IDmat=matrix(rep(1:dim(diffmat)[2],each=dim(diffmat)[1]),nrow=dim(diffmat)[1])
logmat=diffmat>1 & diffmat<threshold
IDfin=IDmat
IDfin[logmat==FALSE]=NA
NegFlux=which(diffmat<threshold^0.2,arr.ind=TRUE)
if(length(NegFlux)>0){
NegFlux[,2]=NegFlux[,2]-1
IDfin[NegFlux]=IDmat[NegFlux]
IDfin[NegFlux[NegFlux[,2]==0,1],1]=0
}
tempout=suppressWarnings(apply(IDfin,1,min,na.rm=TRUE))
tempout[is.infinite(tempout)]=dim(diffmat)[2]
tempout+1
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.