Get the latest version of ProFound and ProFit:
library(devtools)
install_github('asgr/ProFound')
install_github('ICRAR/ProFit')Set global evaluate (basically TRUE for GitHub version and FALSE for CRAN):
evalglobal=TRUEFirst load the libraries we need:
library(ProFound)## Loading required package: Rfits## Loading required package: magicaxis## Loading required package: Rcpplibrary(ProFit)## Loading required package: FITSiolibrary(Rfits)
library(Rwcs)image = Rfits_read_image(system.file("extdata", 'VIKING/mystery_VIKING_Z.fits', package="ProFound"))
profound = profoundProFound(image, magzero=30, verbose=TRUE, plot=TRUE, boundstats=TRUE, rotstats=TRUE)## Running ProFound:## Using raw flux units## Supplied image is 356 x 356 pixels## Extracted pixel scale from header provided: 0.339 asec/pixel## Supplied image is 2.011 x 2.011 amin, 0.001124 deg-sq## Making initial sky map - 0.001 sec## - Sky statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.7399 -0.2272 0.4102 0.5812 1.0496 4.9236## - Sky-RMS statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.159 8.857 9.130 9.134 9.307 10.300## Making initial segmentation image - 0.031 sec## - Running MakeSegim:## - Skipping making initial local estimate of the sky - User provided sky## - Skipping making initial local estimate of the sky RMS - User provided sky RMS## - Smoothing the image - 0.002 sec## - Watershed de-blending - 0.078 sec## - Skipping segmentation plot - plot set to FALSE## - Skipping making final local estimate of the sky - User provided sky## - Skipping making final local estimate of the sky RMS - User provided sky RMS## - Skipping segmentation statistics - segstats set to FALSE or no segments## - MakeSegim is finished! - 0.086 sec## Doing initial aggressive dilation - 0.12 sec## - Running MakeSegimDilate:## - Dilating segments - 0.001 sec## - Skipping segmentation statistics - segstats set to FALSE## - profoundMakeSegimDilate is finished! - 0.016 sec## Making better sky map - 0.139 sec## - Sky statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.1923 -0.7818 -0.1067 -0.1494 0.4264 3.5583## - Sky-RMS statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.246 8.823 8.863 8.903 9.081 9.274## Calculating initial segstats - 0.193 sec## Doing dilations:## Iteration 1 of 6 - 0.211 sec## Iteration 2 of 6 - 0.22 sec## Iteration 3 of 6 - 0.228 sec## Iteration 4 of 6 - 0.236 sec## Iteration 5 of 6 - 0.245 sec## Iteration 6 of 6 - 0.252 sec## - Running MakeSegimDilate:## - Dilating segments - 0 sec## - Skipping segmentation statistics - segstats set to FALSE## - profoundMakeSegimDilate is finished! - 0.033 sec## Doing final aggressive dilation - 0.292 sec## Making final sky map - 0.292 sec## Making final sky RMS map - 0.316 sec## - Sky statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.3946 -0.8394 -0.3054 -0.2772 0.3861 2.8269## - Sky-RMS statistics :## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.247 8.768 8.847 8.871 9.034 9.198## Calculating final segstats for 50 objects - 0.348 sec## - magzero = 30## - gain = NULL (ignored)## - pixscale = 0.339## - rotstats = TRUE## - boundstats = TRUE## - groupstats = TRUE - 0.725 sec## - deblend = FALSE## Plotting segments - 1.102 sec## ProFound is finished! - 1.332 sec
We can identify the large group of central potentially confused sources from the group output from ProFound:
profound$group$groupsegID[1,]## groupID segID Ngroup Npix
## 1 1 1, 2, 3, 6 4 5777The main group has 8 segments, which we can view easily:
magimage(profound$group$groupim==1)
We can check the properties of the group 1 segments too:
group1segstats=profound$segstats[profound$segstats$segID %in% unlist(profound$group$groupsegID[1,"segID"]),]
group1segstats## segID uniqueID xcen ycen xmax ymax RAcen Deccen RAmax
## 1 1 222149 220.8657 148.4511 221.5 148.5 352.2866 -31.82508 352.2866
## 2 2 227142 227.1816 140.4306 226.5 141.5 352.2859 -31.82584 352.2860
## 3 3 194148 193.5578 148.0789 193.5 147.5 352.2897 -31.82512 352.2897
## 6 6 178179 177.3096 178.8072 177.5 178.5 352.2915 -31.82222 352.2914
## Decmax sep flux mag flux_app mag_app cenfrac
## 1 -31.82508 0.2156709 2969722.82 13.81821 571718.03 15.60705 0.031670102
## 2 -31.82574 0.4299177 1004804.67 14.99480 206110.83 16.71475 0.033439675
## 3 -31.82517 0.1972318 56603.27 18.11790 10231.11 19.97519 0.030576097
## 6 -31.82225 0.1225217 45409.81 18.35713 1523.23 22.04309 0.004986304
## N50 N90 N100 R50 R90 R100 SB_N50 SB_N90
## 1 28.74951 146.0390 1858 1.101001 2.481458 8.851060 15.86836 16.99477
## 2 24.17834 141.6071 2061 1.063104 2.572794 9.815253 16.85694 18.13790
## 3 32.73981 202.8898 386 1.185293 2.950651 4.069878 20.30916 21.65144
## 6 184.71767 666.3421 1472 4.480911 8.510611 12.649285 22.42697 23.18176
## SB_N100 xsd ysd covxy corxy con asymm
## 1 19.64182 3.509788 3.543925 1.7629331 0.14173280 0.4436912 0.2338434
## 2 20.93099 4.047922 3.867643 3.7193626 0.23756910 0.4132099 0.3962359
## 3 22.23536 3.209266 3.732403 0.6567635 0.05482956 0.4017057 0.2598715
## 6 23.92789 9.448844 9.427492 -71.0208652 -0.79728042 0.5265087 0.1504434
## flux_reflect mag_reflect semimaj semimin axrat ang signif
## 1 3867533.06 13.53141 3.780130 3.279528 0.8675702 136.95342 Inf
## 2 1890436.91 14.30859 4.420711 3.459526 0.7825723 129.57029 Inf
## 3 72316.08 17.85191 3.758896 3.204309 0.8524601 170.05724 Inf
## 6 44126.83 18.38824 12.656367 4.259262 0.3365312 45.08129 Inf
## FPlim flux_err mag_err flux_err_sky flux_err_skyRMS flux_err_shot
## 1 2.179146 1084.3188 0.0003964290 1012.1471 388.9803 0
## 2 2.137907 781.7274 0.0008446913 666.8667 407.9050 0
## 3 2.742820 213.3895 0.0040931327 122.5519 174.6886 0
## 6 2.269641 1043.6748 0.0249539790 986.3675 341.0810 0
## flux_err_cor cor_seg sky_mean sky_sum skyRMS_mean Nedge Nsky Nobject
## 1 0 0 1.9713314 3662.7338 9.023775 201 49 152
## 2 0 0 2.1602563 4452.2883 8.984926 168 98 70
## 3 0 0 0.9480495 365.9471 8.891345 71 3 68
## 6 0 0 -0.4599773 -677.0866 8.889865 135 107 28
## Nborder Nmask edge_frac edge_excess flag_border iter origfrac Norig
## 1 0 0 0.24378109 1.3104724 0 4 0.9980130 1235
## 2 0 0 0.58333333 1.0322778 0 5 0.9940484 1002
## 3 0 0 0.04225352 1.0145874 0 0 1.0000000 386
## 6 0 0 0.79259259 0.8108876 0 3 0.8864669 637
## skyseg_mean flag_keep
## 1 NA TRUE
## 2 NA TRUE
## 3 NA TRUE
## 6 NA TRUEIn our Full Monty vignette example we detail how to extract a reasonable PSF. We use the solution here.
psf_modellist=list(
moffat=list(
xcen=75/2,
ycen=75/2,
mag=0,
fwhm=3.8,
con=2.04,
ang=0,
axrat=1,
box=0
)
)
psf_model=profitMakeModel(modellist=psf_modellist, dim=c(75,75))$zWe can see where this FWHM lies in comparison to the values estimated from ProFound.
maghist(profound$segstats$R50*2/0.339,breaks=20)## Summary of used sample:## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.617 5.194 6.156 6.925 6.976 26.436## Pop Std Dev: 3.6071## MAD: 1.2541## Half 16-84 Quan (1s): 1.4176## Half 02-98 Quan (2s): 6.056## Using 50 out of 50abline(v=4.47, lty=2)
Have a quick check:
magimage(psf_model)
Let’s try fitting all the sources as Sersic profiles!
group1_modellist=list(
pointsource=list(
xcen= group1segstats$xmax[1:3],
ycen= group1segstats$ymax[1:3],
mag= group1segstats$mag[1:3]
),
sersic=list(
xcen= group1segstats$xmax[4],
ycen= group1segstats$ymax[4],
mag= group1segstats$mag[4],
re= group1segstats$R50[4]/0.339,
nser= 1,
ang= group1segstats$ang[4],
axrat= group1segstats$axrat[4],
box= 0
)
)
group1_interval=list(
pointsource=list(
xcen= rep(list(c(-1,1)),3),
ycen= rep(list(c(-1,1)),3),
mag= rep(list(c(-2,2)),3)
),
sersic=list(
xcen= list(c(-1,1)),
ycen= list(c(-1,1)),
mag= list(c(-2,2)),
re= list(c(0.5,2)),
nser= list(c(0.5,10)),
ang= list(c(-180,360)),
axrat= list(c(0.1,1)),
box= list(c(-1,1))
)
)
for(i in 1:3){
group1_interval$pointsource$xcen[[i]]=group1_modellist$pointsource$xcen[i]+group1_interval$pointsource$xcen[[i]]
group1_interval$pointsource$ycen[[i]]=group1_modellist$pointsource$ycen[i]+group1_interval$pointsource$ycen[[i]]
group1_interval$pointsource$mag[[i]]=group1_modellist$pointsource$mag[i]+group1_interval$pointsource$mag[[i]]
}
for(i in 1){
group1_interval$sersic$xcen[[i]]=group1_modellist$sersic$xcen[i]+group1_interval$sersic$xcen[[i]]
group1_interval$sersic$ycen[[i]]=group1_modellist$sersic$ycen[i]+group1_interval$sersic$ycen[[i]]
group1_interval$sersic$mag[[i]]=group1_modellist$sersic$mag[i]+group1_interval$sersic$mag[[i]]
group1_interval$sersic$re[[i]]=group1_modellist$sersic$re[i]*group1_interval$sersic$re[[i]]
}
group1_tofit=list(
pointsource=list(
xcen= rep(TRUE, 3),
ycen= rep(TRUE, 3),
mag= rep(TRUE, 3)
),
sersic=list(
xcen= TRUE,
ycen= TRUE,
mag= TRUE,
re= TRUE,
nser= TRUE,
ang= TRUE,
axrat= TRUE,
box= FALSE
)
)
group1_tolog=list(
pointsource=list(
xcen= rep(FALSE, 3),
ycen= rep(FALSE, 3),
mag= rep(FALSE, 3)
),
sersic=list(
xcen= FALSE,
ycen= FALSE,
mag= FALSE,
re= TRUE,
nser= TRUE,
ang= FALSE,
axrat= TRUE,
box= FALSE
)
)sigma=profoundMakeSigma(image$imDat, objects=profound$objects, gain=0.5, sky=profound$sky, skyRMS=profound$skyRMS, plot=TRUE)
group1_Data=profitSetupData(image$imDat-profound$sky, sigma=sigma, modellist=group1_modellist, tofit=group1_tofit, tolog=group1_tolog, intervals=group1_interval, magzero=30, algo.func='optim', psf=psf_model, region=matrix(profound$segim %in% unlist(profound$group$groupsegID[1,'segID'])[1:5], 356), verbose=FALSE)profitLikeModel(parm=group1_Data$init, Data=group1_Data, makeplots=TRUE, plotchisq=TRUE)
## [1] -40902.85Next we will run a simple optimisation scheme (because we do not have all day…)
group1_fit=optim(group1_Data$init, profitLikeModel, method='BFGS', Data=group1_Data, control=list(fnscale=-1))And check the results!
profitLikeModel(parm=group1_fit$par, Data=group1_Data, makeplots=TRUE, plotchisq=TRUE)
## [1] -11208.32That is looking pretty pretty pretty good. A bit of residual structure in the central part, but that is as we might expect seeing as we did not attempt to create a more complex (e.g. elliptical and boxy) PSF.
We can compare the inputs and outputs easily too:
group1_modellist_fit=profitRemakeModellist(parm=group1_fit$par, Data=group1_Data)$modellist
group1_modellist_fit## $pointsource
## $pointsource$xcen
## [1] 221.1036 226.3335 193.2357
##
## $pointsource$ycen
## [1] 148.3058 141.0860 147.7205
##
## $pointsource$mag
## [1] 13.82078 14.98410 18.20252
##
##
## $sersic
## $sersic$xcen
## [1] 176.5
##
## $sersic$ycen
## [1] 179.5
##
## $sersic$mag
## [1] 18.40362
##
## $sersic$re
## [1] 13.24699
##
## $sersic$nser
## [1] 0.7775235
##
## $sersic$ang
## [1] 44.04531
##
## $sersic$axrat
## [1] 0.144188
##
## $sersic$box
## [1] 0We can also see the change in a nicely formatted way if we use the relist function in R:
relist(unlist(group1_modellist)-unlist(group1_modellist_fit), skeleton = group1_modellist)The stars have migrated position a small amount, as have the extended sources. The extended sources have both hit a position limit, so it might be sensible to expand the allowed x and y centre intervals. Other parameters are barely changed, e.g. the initial ProFound guesses for the Re (re) and angle (ang) are actually very good. The axial ratios for the extended sources has changed a bit, but this should be expected since these are highly elliptical sources and the convolution with the PSF will be puffing them out when measured without accounting for the PSF in ProFound.
The two bright PSF magnitudes have changed the most, which is not surprising since they were blended together. It is notable that the raw magnitudes calculated by ProFound are actually nearer to the fit values than the rotstats variants. In practice, there is rarely an occasion where the raw segmented magnitudes are not a good enough a starting point for optimisation in ProFit, since it is in practice hard for them to be more than a factor two wrong in flux, which is easily accurate enough for a sensible starting point with ProFit.