ProFit: Making a Complex Model Image
In ProFit: Fit Projected 2D Profiles to Galaxy Images

Get the latest version of ProFit:

library(devtools)
install_github('ICRAR/ProFit')

First load the libraries we need:

library(ProFit)

A Simple Model Structure

It is possible to use the fast low level interface to the profitMakeModel function directly. The example below will generate a 100x100 image matrix, and will run in ~5 microseconds:

ExampleImage0=profitMakeModel(list(sersic=list(xcen=50, ycen=50, re=4, nser=4, ang=30, axrat=0.5)))
str(ExampleImage0)

The structure is a pure image matrix. It is possible to use base image to view this matrix, but the defaults will make it look pretty dreadful, hard to comprehend, and will be slow to actually plot:

image(ExampleImage0)

We can make a much better image of this model easily using magimage, which chooses useful scaling for astronomical images:

magimage(ExampleImage0)

We are not limited to simpel single Sersic profiles, highly complex multi-component systems can also be constructed.

A Multi-Component Structure

First we need to setup a legal model structure, with each component type being a separate list. Whilst this seems more fiddly than before, it allows us to make very complex structures and images very easily (i.e. it scales with complexity very well).

The top level list can be called whatever you like - model is recommended. The model can contain any number of the following components: sersic / coresersic / moffat / king / brokenexp / pointsource / sky. Here is a legal model structure containing 2 Sersic components, 3 stars (point sources) and a sky pedestal:

model1=list(
  sersic=list(
        xcen=c(180, 60),
        ycen=c(90, 10),
        mag=c(15, 13),
        re=c(14, 5),
        nser=c(3, 10),
        ang=c(46, 80),
        axrat=c(0.4, 0.6),
        box=c(0.5,-0.5)
    ),
    pointsource=list(
        xcen=c(34.3,10.1,150),
        ycen=c(74.2,120.5,130.4),
        mag=c(10,13,16)
    ),
    sky=list(
        bg=3e-12
    )
)

With this structure we can make a simple model using profitMakeModel:

ExampleImage1=profitMakeModel(modellist=model1, dim=c(200,200))
str(ExampleImage1)

Next we can make a greyscale image of our requested model structure:

magimage(ExampleImage1)

Notice there are no stars and the models look very sharp - this is because we have not specified a PSF with which to convolve the image. We can use profitPointSource to create a default circular Gaussian PSF with FWHM=2 pixels:

ExamplePSF=profitMakePointSource()
ExamplePSF=ExamplePSF/sum(ExamplePSF)
ExampleImagePSF1=profitMakeModel(modellist=model1, psf=ExamplePSF, dim=c(200,200))

magimage(ExampleImagePSF1)

That looks much more realistic. The convolution adds a modest ~15 ms to the run time, as the PSF is quite small.

You can also specify the PSF as a Sersic model - specifically a Gaussian, where nser=0.5 and re=FWHM/2~1.17sigma:

model1a=model1
model1a$psf=list(sersic=list(nser=0.5,mag=0,re=1,axrat=1,ang=0))
ExampleImagePSF1A=profitMakeModel(modellist=model1a, psf=ExamplePSF, dim=c(200,200))

Note we are still using the PSF image to convolve the Sersic models, but we use the PSF model to create more accurate images of the point sources. This is because we are integrating the flux in the point source more accurately than by interpolating a static PSF image.

We can now easily make a classic bulge+disk system where the two Sersic model components sit on top of each other:

modelBD=list(
    sersic=list(
        xcen=c(100, 100),
        ycen=c(100, 100),
        mag=c(14, 12),
        re=c(2, 15),
        nser=c(4, 1),
        ang=c(0, 60),
        axrat=c(1, 0.3),
        box=c(0,0)
    )
)
ExampleImage=profitMakeModel(modellist=modelBD, psf=ExamplePSF, dim=c(200,200))

magimage(ExampleImage)

The cuspy bulge is evident in the centre of the image, and you can see the dull glow of its large radii wings too (without a truncation high-Nser systems both have a cuspy core and Lorentzian wings that extend off to large radii).

Note that the accuracy of convolution is limited by the pixel dimensions relative to the PSF. In this case, the PSF is quite small, so for accurate convolution, we need to sample the model more finely:

Finesample=3L
FinePSFmodel=model1a$psf
FinePSFmodel$sersic$re=FinePSFmodel$sersic$re*Finesample
ExamplePSFFine=profitMakePointSource(modellist=FinePSFmodel, image=matrix(0,27,27))
ExampleImageFine=profitMakeModel(modellist=modelBD, dim=c(200,200), finesample=Finesample, psf=ExamplePSFFine)

We can see how much of a difference this makes and where:

magimage(ExampleImageFine$z-ExampleImage$z)
maghist((ExampleImageFine$z-ExampleImage$z)/ExampleImageFine$z,xlim=5)

This fine sampling step required manual resampling of the PSF by rescaling re and enlarging the PSF image. When fitting models, profitSetupData will handle these steps automatically.

Make a complex mock image

We can use profitMakeModel to make more complex image structure very easily. The example model below will randomly produce 20 Sersic profiles, 10 PSFs and a sky background:

model2=list(
    sersic=list(
        xcen=runif(20,0,200),
        ycen=runif(20,0,200),
        mag=runif(20,15,20),
        re=runif(20,1,100),
        nser=runif(20,0.5,8),
        ang=runif(20,0,180),
        axrat=runif(20,0.3,1),
        box=runif(20,-0.3,0.3)
    ),
    pointsource=list(
        xcen=runif(10,0,200),
        ycen=runif(10,0,200),
        mag=runif(10,15,20)
    ),
    sky=list(
        bg=3e-12
    )
)

As before we can run this through profitMakeModel and then plot the image. This should run in ~0.5 second, which means we are scaling well with our more complex model:

system.time(ExampleImagePSF2<-profitMakeModel(modellist=model2, psf=ExamplePSF, dim=c(200,200)))
magimage(ExampleImagePSF2, hi=1)

We can try a bigger image (more like a typical astronomy image) but this will take longer to generate (~3 seconds). However, this is much faster than the naive 1x(1000/200)^2=25 seconds that we might expect. This is because the lower level C++ routine adaptively adjusts to calculate the pixel integrals. Pixels near the centre need more subdivisions to accurately determine the flux, whereas at large radii typically fewer calculations are required. Hence generating a much bigger image does not directly scale the computation time.

model3=list(
    sersic=list(
        xcen=runif(20,0,1000),
        ycen=runif(20,0,1000),
        mag=runif(20,15,20),
        re=runif(20,1,100),
        nser=runif(20,0.5,8),
        ang=runif(20,0,180),
        axrat=runif(20,0.3,1),
        box=runif(20,-0.3,0.3)
    ),
    pointsource=list(
        xcen=runif(10,0,1000),
        ycen=runif(10,0,1000),
        mag=runif(10,15,20)
    ),
    sky=list(
        bg=3e-12
    )
)
system.time(ExampleImagePSF3<-profitMakeModel(modellist=model3, psf=ExamplePSF, dim=c(1000,1000)))
magimage(ExampleImagePSF3)

A faster image can be generated if the remax flag is turned on. With remax=10 it only takes ~2 seconds to make the model image. This might be accurate enough for many simulation purposes and can be a factor of 2 faster.

system.time(profitMakeModel(modellist=model3, psf=ExamplePSF, dim=c(1000,1000), remax=10))

It is easy to add together model images using `profitAddMats', giving flexibility in how mock images are made:

ExampleImageAdd=profitAddMats(ExampleImagePSF3$z, ExampleImagePSF2$z, c(300,400))
magimage(ExampleImageAdd)

In this example it is evident that there is a lot of low surface brightness flux in `ExampleImagePSF2', so a much bigger image should be made to remove the flux discontinuity at the edges of the 200x200 sub-image.