ProFound: It’s All in the Blend

Aaron Robotham

2025-11-12

evalglobal=TRUE
library(ProFound)
#> Loading required package: Rfits
#> Loading required package: magicaxis
#> Loading required package: Rcpp
#> 
#> Attaching package: 'ProFound'
#> The following object is masked _by_ '.GlobalEnv':
#> 
#>     %nin%
library(ProFit)
library(RANN)
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

Over the years ProFound has developed a few distinct and competing approaches to source de-blending. The core profoundProFound function contains a few different options to achieve this (see the deblendtype argument options), but here we will use the “fit” mode that attempts to fit a simple profile to each blended object and use the approximate profile as a de-blending kernel. This has been present in ProFoun for many years, but is rarely used in practice (perhaps people just do not know the option is there?) In practice, as we shall see, it is often not much better than the watershed based de-blending scheme we use in the base function, even though the latter creates hard de-blending boundaries and does not share the flux of a pixel between segments. This seems surprising, but is a natural consequence of what the saddle point in flux usually represents.

In late 2025 I developed a new elliptical aperture routine built in Rcpp that is extremely fast (faster than the equivalent SExtractor code) and can even be threaded with OpenMP. As well as offering proper circular and elliptical aperture photometry, it offers a few useful approaches for de-blending sources across hard pixel bouandaries, i.e. it can fractional distribute flux between an arbitrarily large number of sources that appear to overlap. Much like SExtractor, we can estimate the properties of the ellipses to use from the naive segment-based profoundProFound source extraction.

Resolved Data Example

First let us make. a mock image where we have a ground truth for some reasonable looking galaxies. To do this we will use ProFit, one of the core **ProTools* software packages that is used for generating and fitting 2D light profiles of many sources.

ExamplePSF = profitMakeGaussianPSF()

set.seed(1)
model = list(
  sersic = list(
    xcen = runif(50, 100, 900),
    ycen = runif(50, 100, 900),
    mag = runif(50, 15, 20),
    re = runif(50, 3, 10),
    nser = runif(50, 1, 4),
    ang = runif(50, 0, 180),
    axrat = runif(50, 0.3, 1),
    box = runif(50, -0.3, 0.3)
  ),
  pointsource = list(
    xcen = runif(20, 100, 900),
    ycen = runif(20, 100, 900),
    mag = runif(20, 15, 20)
  )
)

model_core = data.frame(xcen = c(model$sersic$xcen, model$pointsource$xcen),
                        ycen = c(model$sersic$ycen, model$pointsource$ycen),
                        mag = c(model$sersic$mag, model$pointsource$mag)
                        )

ExampleImage = profitMakeModel(modellist=model, psf=ExamplePSF, dim=c(1000,1000))$z

ExampleImage = ExampleImage + rnorm(1e6, sd=1e-12)

First let us run a vanilla version of profoundProFound on this and see how we do:

pro_out = profoundProFound(ExampleImage, deblend=TRUE, plot=TRUE)

Find the matching sources:

pro_match = nn2(pro_out$segstats[,c("xcen", "ycen")], model_core[,c("xcen","ycen")],
                searchtype='radius', radius=2)
magplot(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag"],
        ylim=c(-0.2,0.2), xlab='True [mag]', ylab='ProFound - True [mag]')
points(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag_db"],
       col='red')
legend('topleft', legend=c('Segment','De-blend'), col=c('black','red'), pch=1)

Now we will try the newer elliptical based routine in ProFound. Notice we still need various outputs from the standard run of profoundProFound to estimate a good de-blend solution. We could not supply the rad_re and the wt, but the solution would be much less good (try this yourself!).

First, for fun, we ccan look at the weight map that is used internally with and without profile scaling:

ellip_ref = pro_out$segstats[,c("flux", "xcen", "ycen", "R50", "R100",  "ang", "axrat")]

magimage(ExampleImage) #for reference to see where profoundProFound is putting the apertures


magimage(profoundEllipWeight(ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100, ellip_ref$ang, ellip_ref$axrat, dimx = 1000, dimy = 1000))


magimage(profoundEllipWeight(ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100, ellip_ref$ang, ellip_ref$axrat, dimx = 1000, dimy = 1000, wt=ellip_ref$flux, rad_re=ellip_ref$R50))

ellip_flux_deblend =  profoundEllipFlux(ExampleImage, ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100, ellip_ref$ang, ellip_ref$axrat, deblend=TRUE, wt=ellip_ref$flux, rad_re=ellip_ref$R50)
ellip_mag_deblend = profoundFlux2Mag(ellip_flux_deblend)

Re the above routine, it is a very focussed Rcpp level function (basically all C code) that is highly optimised to do elliptical aperture photometry and de-blending. The original profoundProFound routine takes about 2-4 seconds to run even on multiple cores, whereas on a single core profoundEllipFlux takes about 0.1 seconds. It can also be run with full OpenMp threading if the ProFound package is built with OpenMP support (that is beyond the scope of this vignette, but there is a lot of information online about how to build Rcpp functions and packges with OpenMP threading).

Now we can compare all three approaches to the ground truth:

magplot(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag"],
        ylim=c(-0.2,0.2), xlab='True [mag]', ylab='True - ProFound [mag]')
points(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag_db"],
       col='red')
points(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend[pro_match$nn.idx[,1]],
       col='blue')
legend('topleft', legend=c('Segment','Pro De-blend', 'Ellip De-blend'), col=c('black','red','blue'), pch=1)

Perhaps it is better to look at the 1D mag difference density directly:

magplot(density(model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag"]),
        xlab='True - ProFound[mag]', ylab='PDF', xlim=c(-0.2,0.2), ylim=c(0,23))
lines(density(model_core[pro_match$nn.idx[,1] > 0,"mag"] - pro_out$segstats[pro_match$nn.idx[,1],"mag_db"]), col='red')
lines(density(model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend[pro_match$nn.idx[,1]]), col='blue')
legend('topleft', legend=c('Segment','Pro De-blend', 'Ellip De-blend'), col=c('black','red','blue'), pch=1)

There is not a lot in it, but I would argue based on the above the newer Elliptical de-blending scheme does a little better than either the standard profoundProFound segment based scheme, or the ‘fit’ method of de-blending that is available within the function. So I guess that’s the conclusion- use it if you want better de-blending.

It’s interesting to note that expanding the apertures has little affect on the flux with the de-blending being used in this way, since so little distant flux will be claimed by an expanded aperture. This suggests that when used in this mode, users should consider inflating their apertures a bit to guarantee all the available flux is attributed to each object. The caveat there is you start adding noise, but for extended sources the bigger error tends to be aperture systematics rather than sky noise.

As you can see, the nominal apertures now massively overlap on sky, but the de-blending weighting helps to control this:

magimage(ExampleImage) #for reference to see where profoundProFound is putting the apertures


magimage(profoundEllipWeight(ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100*2, ellip_ref$ang, ellip_ref$axrat, dimx = 1000, dimy = 1000))


magimage(profoundEllipWeight(ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100*2, ellip_ref$ang, ellip_ref$axrat, dimx = 1000, dimy = 1000, wt=ellip_ref$flux, rad_re=ellip_ref$R50))

ellip_flux_deblend_big =  profoundEllipFlux(ExampleImage, ellip_ref$xcen, ellip_ref$ycen,
  ellip_ref$R100*2, ellip_ref$ang, ellip_ref$axrat, deblend=TRUE, wt=ellip_ref$flux, rad_re=ellip_ref$R50)
ellip_mag_deblend_big = profoundFlux2Mag(ellip_flux_deblend_big)

Overall the performance of the larger apertures is very similar:

magplot(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend[pro_match$nn.idx[,1]],
        xlim=c(15,21), ylim=c(-0.2,0.2), col='blue', xlab='ProFound - True [mag]', ylab='PDF')
points(model_core[pro_match$nn.idx[,1] > 0,"mag"], model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend_big[pro_match$nn.idx[,1]],
       col='darkgreen')
legend('topleft', legend=c('Ellip De-blend', 'Ellip De-blend Big'), col=c('blue','darkgreen'), pch=1)

magplot(density(model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend[pro_match$nn.idx[,1]]),
        col='blue', xlab='ProFound - True [mag]', ylab='PDF', xlim=c(-0.2,0.2), ylim=c(0,23))
lines(density(model_core[pro_match$nn.idx[,1] > 0,"mag"] - ellip_mag_deblend_big[pro_match$nn.idx[,1]]), col='darkgreen')
legend('topleft', legend=c('Ellip De-blend', 'Ellip De-blend Big'), col=c('blue','darkgreen'), pch=1)

Unresolved Data Examples

FitMagPSF Example Data

When fitting totally unresolved data, ProFound already comes with profoundFitMagPSF. This correctly deblends flux into objects using the known PSF of the image.

s250_im = Rfits_read_image(system.file("extdata", 'IRdata/s250_im.fits',
                                       package="ProFound"))
s250_psf = Rfits_read_image(system.file("extdata",'IRdata/s250_psf.fits',
                                        package="ProFound"))$imDat
s250_psf[s250_psf < 0.001] = 0
s250_psf = s250_psf[41:61, 41:61]

magzero_s250 = 11.68

pro_s250 = profoundProFound(s250_im, pixcut=1, skycut=2.5, ext=1, redosky=FALSE, iters=1,
                          tolerance=0, sigma=0, magzero=magzero_s250)
#> dim(image)[1]/box[1] must be >=3, box[1] modified to 40
#> dim(image)[2]/box[2] must be >=3, box[2] modified to 40
#> dim(image)[1]/grid[1] must be >=3, grid[1] modified to 40
#> dim(image)[2]/grid[2] must be >=3, grid[2] modified to 40
pro_s250$segstats = pro_s250$segstats[!is.na(pro_s250$segstats$mag),]

newmag = profoundFitMagPSF(RAcen=pro_s250$segstats$RAcen, Deccen=pro_s250$segstats$Deccen,
                         image=s250_im, psf=s250_psf, doProFound=TRUE, findextra=TRUE,
                         verbose=TRUE, redosky=FALSE, magzero=magzero_s250)
#> Running initial ProFound
#> dim(image)[1]/box[1] must be >=3, box[1] modified to 40
#> dim(image)[2]/box[2] must be >=3, box[2] modified to 40
#> dim(image)[1]/grid[1] must be >=3, grid[1] modified to 40
#> dim(image)[2]/grid[2] must be >=3, grid[2] modified to 40
#> - using ProFound skyRMS model for im_sigma
#> Iterating over source model
#> - iteration 1 of 5
#> - iteration 2 of 5
#> - iteration 3 of 5
#> - iteration 4 of 5
#> - iteration 5 of 5
#> Finding extra sources with ProFound
#> dim(image)[1]/box[1] must be >=3, box[1] modified to 40
#> dim(image)[2]/box[2] must be >=3, box[2] modified to 40
#> dim(image)[1]/grid[1] must be >=3, grid[1] modified to 40
#> dim(image)[2]/grid[2] must be >=3, grid[2] modified to 40
#> - gained 2 additional sources
#> Rerunning source model with extra sources
#> Iterating over source model
#> - iteration 1 of 5
#> - iteration 2 of 5
#> - iteration 3 of 5
#> - iteration 4 of 5
#> - iteration 5 of 5
magimage(s250_im$imDat - newmag$finalmodel, qdiff=T)

It is worth noting how similar a job we can do using the new profoundAperFlux function that makes use of profile based source de-blending.

First we need to specifying some properties of the PSF:

R50 = 1.7 #Re in pixels
R100 = 5 #Maximum radius to use for the apertures

magimage(s250_im$imDat, qdiff=TRUE)


magimage(profoundAperWeight(newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, dimx = 120, dimy = 120))


magimage(profoundAperWeight(newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, dimx = 120, dimy = 120, rad_re=R50, nser=0.5))

sky = sum(s250_im$imDat - newmag$finalmodel)/(120^2)

newflux_db = profoundAperFlux(s250_im$imDat - sky, newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, deblend=TRUE)
newmag_db = profoundFlux2Mag(newflux_db, magzero=magzero_s250)

newflux_re = profoundAperFlux(s250_im$imDat - sky, newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, rad_re=R50, nser=0.5, deblend=TRUE)
newmag_re = profoundFlux2Mag(newflux_re, magzero=magzero_s250)

wt = newmag$psfstats$flux
wt[is.na(wt)] = 0
newflux_wt = profoundAperFlux(s250_im$imDat - sky, newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, rad_re=R50, nser=0.5, wt=wt, deblend=TRUE)
newmag_wt = profoundFlux2Mag(newflux_wt, magzero=magzero_s250)

And now we can compare how it looks:

magplot(newmag$psfstats$mag, newmag$psfstats$mag - newmag_wt, xlim=c(11,17), ylim=c(-2,2),
        xlab='FitMagPSF [mag]', ylab='PSF - DB [mag]')

We can even make a reasonable difference image using the weight map:

magimage(s250_im$imDat, qdiff=TRUE)


magimage(profoundAperWeight(newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, dimx = 120, dimy = 120, rad_re=R50, nser=0.5, wt=newflux_wt))


magimage(s250_im$imDat - profoundAperWeight(newmag$psfstats$xcen, newmag$psfstats$ycen,
  R100, dimx = 120, dimy = 120, rad_re=R50, nser=0.5, wt=newflux_wt/3), qdiff=TRUE)

Where the new profoundAperFlux function could be advantageous is when you have a mixture of resolved and unresolved sources in your image. The pure PSF modelling profoundFitMagPSF cannot properly account for partially resolved sources, and in general they will be biased to having erroneously low flux.

Mock Data

Now we are going to create a crowded unresolved field that normally we would use profoundFitMagPSF on, but use our new aperture photometry and deblending code, where we how have a ground-truth to compare to.

set.seed(1)
model = list(
  pointsource = list(
    xcen = runif(800, 10, 190),
    ycen = runif(800, 10, 190),
    mag = runif(800, 15, 20)
  )
)

model_core = data.frame(xcen = model$pointsource$xcen,
                        ycen = model$pointsource$ycen,
                        mag = model$pointsource$mag
                        )

ExampleImage = profitMakeModel(modellist=model, psf=s250_psf, dim=c(200,200))$z

ExampleImage = ExampleImage + rnorm(200^2, sd=1e-8)

Let’s take a look:

magimage(ExampleImage)

Great, looks nice and confusing :-D

We are going to assume a best case scenario that somehow we know where the sources are, and we will run profoundFitMagPSF and profoundAperFlux at these positions:

mag_FitMag = profoundFitMagPSF(xcen=model_core$xcen, ycen=model_core$ycen,
               im_sigma=matrix(1e-8, 200, 200), image=ExampleImage, psf=s250_psf,
               verbose=TRUE, redosky=FALSE)
#> Iterating over source model
#> - iteration 1 of 5
#> - iteration 2 of 5
#> - iteration 3 of 5
#> - iteration 4 of 5
#> - iteration 5 of 5

Notice rad_re below is the approximate value for the PSF used in this work (so replace as appropriate). Our starting point for the weight map looks like this:

wtmap = profoundAperWeight(cx=model_core$xcen, cy=model_core$ycen, rad=5, rad_re=1.7,
                           dimx=200, dimy=200)
magimage(wtmap)

flux_AperFlux = profoundAperFlux(ExampleImage, cx=model_core$xcen, cy=model_core$ycen,
                                rad=5, rad_re=1.7, deblend=TRUE, iterations=5)
mag_AperFlux = profoundFlux2Mag(flux_AperFlux)

How do they compare to the ground truth?

magplot(model_core$mag, model_core$mag - mag_FitMag$psfstats$mag, asp=1, xlim=c(15,20),
        ylim=c(-2,2), xlab='True [mag]', ylab='True - Measure [mag]')
points(model_core$mag, model_core$mag - mag_AperFlux, col='red')
legend('topleft', legend=c('FitMagPDF', 'AperFlux'), col=c('black', 'red'), pch=1)

And as density plots:

magplot(density(model_core$mag - mag_FitMag$psfstats$mag, na.rm=TRUE), xlab='True - Measure [mag]', ylab='PDF', xlim=c(-2,2), ylim=c(0,1.3))
lines(density(model_core$mag - mag_AperFlux, na.rm=TRUE), col='red')
legend('topleft', legend=c('FitMagPDF', 'AperFlux'), col=c('black','red'), pch=1)

So a full profoundFitMagPSF produces better photometry (as we should expect), but we get fairly similar performance from profoundAperFlux. It is also much faster to run (about x10), so in certain use cases that might be th deciding factor. profoundAperFlux is also a very low level function (literally Rcpp code) so it lacks some of the friendly features we often expect, and it does not produce native error for the fluxes. For aperture errors you often do best throwing down random apertures in any case, based on the emptiest positions in the profoundAperWeight map.

The example above was a bit artificial in that we pretended we knew the positions of the sources, but that’s not too far from the truth in many cases. We often have very deep NIR data when we also have confused and unresolved FIT data, and the former is a very good starting point for the later.