This is the highest level source detection function provided in ProFit, calculating both the initial segmentation map and reasonable estimates for the total flux apertures for each source in an automatic manner.
profoundProFound(image, segim, objects, mask, skycut=1, pixcut=3, tolerance = 4, ext = 2,
sigma = 1, smooth = TRUE, SBlim, size = 5, shape = "disc", iters = 6, threshold = 1.05,
converge = 'flux', magzero = 0, gain = NULL, pixscale = 1, sky, skyRMS, redosky = TRUE,
redoskysize = 21, box = c(100,100), grid = box, type = "bilinear", skytype = "median",
skyRMStype = "quanlo", sigmasel = 1, doclip = TRUE, shiftloc = FALSE, paddim = TRUE,
header, verbose = FALSE, plot = FALSE, stats = TRUE, rotstats = FALSE, boundstats = FALSE,
nearstats=boundstats, groupstats=boundstats, offset = 1, sortcol = "segID",
decreasing = FALSE, lowmemory=FALSE, keepim = TRUE, ...)Numeric matrix; required, the image we want to analyse. If image is a list as created by readFITS, read.fits of magcutoutWCS then the image part of these lists is passed to image and the correct header part is passed to header. Note, image NAs are treated as masked pixels.
Numeric matrix; a specified segmentation map of the image. This matrix *must* be the same dimensions as image if supplied. If this is option is used then profoundProFound will not compute its initial segmentation map using profoundMakeSegim, which is then dilated. Instead it will use the one passed through segim.
Boolean matrix; optional, object mask where 1 is object and 0 is sky. If provided, this matrix *must* be the same dimensions as image.
Boolean matrix; optional, parts of the image to mask out (i.e. ignore), where 1 means mask out and 0 means use for analysis. If provided, this matrix *must* be the same dimensions as image.
Numeric scalar; the lowest threshold to make on the image in units of the skyRMS. Passed to profoundMakeSegim.
Integer scalar; the number of pixels required to identify an object. Passed to profoundMakeSegim.
Numeric scalar; the minimum height of the object in the units of skyRMS between its highest point (seed) and the point where it contacts another object (checked for every contact pixel). If the height is smaller than the tolerance, the object will be combined with one of its neighbours, which is the highest. The range 1-5 offers decent results usually. Passed to profoundMakeSegim.
Numeric scalar; radius of the neighbourhood in pixels for the detection of neighbouring objects. Higher value smooths out small objects. Passed to profoundMakeSegim.
Numeric scalar; standard deviation of the blur used when smooth=TRUE. Passed to profoundMakeSegim.
Logical; should smoothing be done on the target image? Passed to profoundMakeSegim.
Numeric scalar; the mag/asec^2 surface brightness threshold to apply. This is always used in conjunction with skycut, so set skycut to be very large (e.g. Inf) if you want a pure surface brightness threshold for the segmentation. magzero and pixscale must also be present for this to be used. Passed to profoundMakeSegim.
Integer scalar; the size (e.g. width/diameter) of the dilation kernel in pixels. Should be an odd number else will be rounded up to the nearest odd number. See makeBrush. Passed to profoundMakeSegimDilate.
Character scalar; the shape of the dilation kernel. See makeBrush. Passed to profoundMakeSegimDilate.
Integer scalar; the maximum number of curve of growth dilations should be made. This needs to be large enough to capture all the flux for sources of interest, but increasing this will increase the computation time for profoundProFound. If this is set to zero then the initial segim image wither provided or computed internally via profoundMakeSegim will be used instead.
Numeric scalar; After the curve of growth dilations, threshold is the relative change of the converging property (see converge) that flags convergence. If consecutive iterations have a relative difference within this ratio then the dilation is stopped, and this iteration is used to define the segmentation of the object. The effect of this is that different objects will be dilated for a different number of iterations. Usually fainter sources require more.
Character scalar; the segmentation property to compare for relative convergence. The options are in principle any column that is output by profoundSegimStats, but in practice it should be something that increases slowly with dilation and tends to converge when the total flux is being captured. Good options are therefore 'flux' (default), 'R50' and 'R90'.
Numeric scalar; the magnitude zero point. What this implies depends on the magnitude system being used (e.g. AB or Vega). If provided along with pixscale then the flux and surface brightness outputs will represent magnitudes and mag/asec^2.
Numeric scalar; the gain (in photo-electrons per ADU). This is only used to compute object shot-noise component of the flux error (else this is set to 0).
Numeric scalar; the pixel scale, where pixscale=asec/pix (e.g. 0.4 for SDSS). If set to 1 (default), then the output is in terms of pixels, otherwise it is in arcseconds. If provided along with magzero then the flux and surface brightness outputs will represent magnitudes and mag/asec^2.
User provided estimate of the absolute sky level. If this is not provided then it will be computed internally using profoundMakeSkyGrid. Can be a scalar or a matrix matching the dimensions of image (allows values to vary per pixel). This will be subtracted off the image internally, so only provide this if the sky does need to be subtracted!
User provided estimate of the RMS of the sky. If this is not provided then it will be computed internally using profoundMakeSkyGrid. Can be a scalar or a matrix matching the dimensions of image (allows values to vary per pixel).
Logical; should the sky and sky RMS grids be re-computed using the final segmentation map? This uses profoundMakeSkyGrid to compute the sky and sky RMS grids. If redosky=TRUE then the output will include the aggressively masked objects_redo image, if redosky=FALSE then objects_redo will be NA.
Integer scalar; the size (e.g. width/diameter) of the dilation kernel in pixels to apply to the object mask before performing the initial and final aggressively masked sky estimates (the latter is only relevant if redosky=TRUE). Should be an odd number else will be rounded up to the nearest odd number. See makeBrush. Dilation is done by profoundMakeSegimDilate. If redosky=TRUE, the final dilated objects mask is returned as objects_redo. As a rule of thumb you probably want ~50% of your image pixels to be masked as objects, much more than this and you might not be able to sample enough sky pixels, much more less and the sky estimates might be biased by object flux in the wings.
Integer vector; the dimensions of the box car filter to estimate the sky with.
Integer vector; the resolution of the background grid to estimate the sky with. By default this is set to be the same as the box.
Character scalar; either "bilinear" for bilinear interpolation (default) or "bicubic" for bicubic interpolation. The former is safer, especially near edges where bicubic interpolation can go a bit crazy.
Character scalar; the type of sky level estimator used. Allowed options are 'median' (the default), 'mean' and 'mode' (see profoundSkyEstLoc for an explanation of what these estimators do). In all cases this is the estimator applied to unmasked and non-object pixels. If doclip=TRUE then the pixels will be dynamically sigma clipped before the estimator is run.
Character scalar; the type of sky level estimator used. Allowed options are 'quanlo' (the default), 'quanhi', 'quanboth', and 'sd' (see profoundSkyEstLoc for an explanation of what these estimators do). In all cases this is the estimator applied to unmasked and non-object pixels. If doclip=TRUE then the pixels will be dynamically sigma clipped before the estimator is run.
Numeric scalar; the quantile to use when trying to estimate the true standard-deviation of the sky distribution. If contamination is low then the default of 1 is about optimal in terms of S/N, but you might need to make the value lower when contamination is very high.
Logical; should the unmasked non-object pixels used to estimate to local sky value be further sigma-clipped using magclip? Whether this is used or not is a product of the quality of the objects extraction. If all detectable objects really have been found and the dilated objects mask leaves only apparent sky pixels then an advanced user might be confident enough to set this to FALSE. If an doubt, leave as TRUE.
Logical; should the cutout centre for the sky shift from loc of the desired box size extends beyond the edge of the image? (See magcutout for details).
Logical; should the cutout be padded with image data until it meets the desired box size (if shiftloc is true) or padded with NAs for data outside the image boundary otherwise? (See magcutout for details).
Full FITS header in table or vector format. If this is provided then the segmentations statistics table will gain RAcen and Decen coordinate outputs. Legal table format headers are provided by the read.fitshdr function or the hdr list output of read.fits in the astro package; the hdr output of readFITS in the FITSio package or the header output of magcutoutWCS. Missing header keywords are printed out and other header option arguments are used in these cases. See magWCSxy2radec.
Logical; should verbose output be displayed to the user? Since big image can take a long time to run, you might want to monitor progress.
Logical; should a diagnostic plot be generated? This is useful when you only have a small number of sources (roughly a few hundred). With more than this it can start to take a long time to make the plot!
Logical; should statistics on the segmented objects be returned? If magzero and pixscale have been provided then some of the outputs are computed in terms of magnitude and mag/asec^2 rather than flux and flux/pix^2 (see Value).
Logical; if TRUE then the asymm, flux_reflect and mag_reflect are computed, else they are set to NA. This is because they are very expensive to compute compared to other photometric properties.
Logical; if TRUE then various pixel boundary statistics are computed (Nedge, Nsky, Nobject, Nborder, edge_frac, edge_excess and FlagBorder). If FALSE these return NA instead (saving computation time).
Logical; if TRUE then the IDs of nearby segments is calculated via profoundSegimNear and output to the returned object near. By default this option is linked to boundstats, i.e. it is assumed if you want boundary statistics then you probably also want nearby object IDs returned.
Logical; if TRUE then the IDs of grouped segments is calculated via profoundSegimGroup and output to the returned object group. By default this option is linked to boundstats, i.e. it is assumed if you want boundary statistics then you probably also want grouped object IDs returned.
Integer scalar; the distance to offset when searching for nearby segments (used in both profoundSegimStats and profoundSegimNear).
Character; name of the output column that the returned segmentation statistics data.frame should be sorted by (the default is segID, i.e. segment order). See below for column names and contents.
Logical; if FALSE (default) the segmentation statistics data.frame will be sorted in increasing order, if TRUE the data.frame will be sorted in decreasing order.
Logical; if TRUE then a low memory mode of ProFound will be used. This limits the large image pixel matched outputs to just segim, with segim_orig, objects and objects_redo set to NULL, and sky and skyRMS set to 0. Internally the sky and skyRMS are used as normal for flux estimates, but they are removed as soon as possible within the function in order to free up memory.
Logical; if TRUE then the input image and mask matrices are passed through to the image output of the function. If FALSE then this is set to NULL.
Further arguments to be passed to magimage. Only relevant is plot=TRUE.
A object list of class 'profound' containing:
Integer matrix; the dilated and converged segmentation map matched pixel by pixel to image.
Integer matrix; the pre-dilated segmentation map matched pixel by pixel to image.
Logical matrix; the object map matched pixel by pixel to image. 1 means there is an object at this pixel, 0 means it is a sky pixel. Can be used as a mask in various other functions that require objects to be masked out.
Logical matrix; the dilated object map matched pixel by pixel to image. See redosky and redoskysize. Can be used as a mask in various other functions that require objects to be masked out.
The estimated sky level of the image.
The estimated sky RMS of the image.
The input image matrix if keepim=TRUE, else NULL.
The input mask matrix if keepim=TRUE, else NULL.
If stats=TRUE this is a data.frame (see below), otherwise NULL.
The total number of segments extracted (dim(segstats)[1]).
If nearstats=TRUE then contains the output of profoundSegimNear.
If groupstats=TRUE then contains the output of profoundSegimGroup.
The header provided, if missing this is NULL.
The surface brightness limit of detected objects. Requires at least magzero to be provided and skycut>0, else NULL. profoundMakeSegimExpand only.
The assumed magnitude zero point. This is relevant to various outputs returned by the segmentation statistics.
The dimensions of the processed image.
The assumed pixel scale. This is relevant to various outputs returned by the segmentation statistics.
The assumed image gain (if NULL it was not used). This is relevant to various outputs returned by the segmentation statistics.
The original function call.
Segmentation ID, which can be matched against values in segim
Unique ID, which is fairly static and based on the xmax and ymax position
Flux weighted x centre
Flux weighted y centre
x position of maximum flux
y position of maximum flux
Flux weighted degrees Right Ascension centre (only present if a header is provided)
Flux weighted degrees Declination centre (only present if a header is provided)
Right Ascension of maximum flux (only present if a header is provided)
Declination of maximum flux (only present if a header is provided)
Radial offset between the cen and max definition of the centre (units of pixscale, so if pixscale represents the standard asec/pix this will be asec)
Total flux (calculated using image-sky) in ADUs
Total flux converted to mag using magzero
Fraction of flux in the brightest pixel
Number of brightest pixels containing 50% of the flux
Number of brightest pixels containing 90% of the flux
Total number of pixels in this segment, i.e. contains 100% of the flux
Approximate elliptical semi-major axis containing 50% of the flux (units of pixscale, so if pixscale represents the standard asec/pix this will be asec)
Approximate elliptical semi-major axis containing 90% of the flux (units of pixscale, so if pixscale represents the standard asec/pix this will be asec)
Approximate elliptical semi-major axis containing 100% of the flux (units of pixscale, so if pixscale represents the standard asec/pix this will be asec)
Mean surface brightness containing brightest 50% of the flux, calculated as flux*0.5/N50 (if pixscale has been set correctly then this column will represent mag/asec^2. Otherwise it will be mag/pix^2)
Mean surface brightness containing brightest 90% of the flux, calculated as flux*0.9/N90 (if pixscale has been set correctly then this column will represent mag/asec^2. Otherwise it will be mag/pix^2)
Mean surface brightness containing all of the flux, calculated as flux/N100 (if pixscale has been set correctly then this column will represent mag/asec^2. Otherwise it will be mag/pix^2)
Weighted standard deviation in x (always in units of pix)
Weighted standard deviation in y (always in units of pix)
Weighted covariance in xy (always in units of pix)
Weighted correlation in xy (always in units of pix)
Concentration, R50/R90
180 degree flux asymmetry (0-1, where 0 is perfect symmetry and 1 complete asymmetry)
Flux corrected for asymmetry by doubling the contribution of flux for asymmetric pixels (defined as no matching segment pixel found when the segment is rotated through 180 degrees)
flux_reflect converted to mag using magzero
Weighted standard deviation along the major axis, i.e. the semi-major first moment, so ~2 times this would be a typical major axis Kron radius (always in units of pix)
Weighted standard deviation along the minor axis, i.e. the semi-minor first moment, so ~2 times this would be a typical minor axis Kron radius (always in units of pix)
Axial ratio as given by min/maj
Orientation of the semi-major axis in degrees. This has the convention that 0= | (vertical), 45= \, 90= - (horizontal), 135= /, 180= | (vertical)
Approximate singificance of the detection using the Chi-Square distribution
Approximate false-positive significance limit below which one such source might appear spuriously on an image this large
Estimated total error in the flux for the segment
Estimated total error in the magnitude for the segment
Sky subtraction component of the flux error
Sky RMS component of the flux error
Object shot-noise component of the flux error (only if gain is provided)
Mean flux of the sky over all segment pixels
Total flux of the sky over all segment pixels
Mean value of the sky RMS over all segment pixels
Number of edge segment pixels that make up the outer edge of the segment
Number of edge segment pixels that are touching sky
Number of edge segment pixels that are touching another object segment
Number of edge segment pixels that are touching the image border
Number of edge segment pixels that are touching a masked pixel (note NAs in image are also treated as masked pixels)
Fraction of edge segment pixels that are touching the sky i.e. NskyNedge, higher generally meaning more robust segmentation statistics
Ratio of the number of edge pixels to the expected number given the elliptical geometry measurements of the segment. If this is larger than 1 then it is a sign that the segment geometry is irregular, and is likely a flag for compromised photometry
A binary flag telling the user which image borders the segment touches. The bottom of the image is flagged 1, left=2, top=4 and right=8. A summed combination of these flags indicate the segment is in a corner touching two borders: bottom-left=3, top-left=6, top-right=12, bottom-right=9.
The iteration number when the source was flagged as having convergent flux
The ratio between the final converged flux and the initial profoundMakeSegim iso-contour estimate
A suggested flag for selecting good objects. Objects flagged FALSE have hit the iteration limit and have grown their flux by more than the median for all objects at the iteration limit.
This high level function is both a source detection and a segmented aperture growing function. The latter is achieved through consecutive dilation and flux measurement operations. It is not super fast, but it is designed to be fairly robust and fast enough for most use cases.
profoundProFound initially makes a segmentation map using the profoundMakeSegim function. It then makes repeated dilations and flux measurements of this segmentation map using profoundMakeSegimDilate, and calculates the convergent flux segment for each source. These are combined to make a final segmentation map with associated source statistics (if requested).
The defaults should work reasonably well on modern survey data (see Examples), but should the solution not be ideal try modifying these parameters (in order of impact priority): skycut, pixcut, tolerance, sigma, ext.
profoundMakeSegimDilate is similar in nature to the pixel growing objmask routine in IRAF (see the ngrow and agrow description at http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?objmasks). This similarity was discovered after implementation, but it is worth noting that the higher level curve of growth function profoundProFound is not trivially replicated by other astronomy tools.
profoundMakeSegim, profoundMakeSegimDilate, profoundMakeSegimExpand, profoundMakeSegimPropagate, profoundSegimStats, profoundSegimPlot
# NOT RUN {
image=readFITS(system.file("extdata", 'VIKING/mystery_VIKING_Z.fits', package="ProFound"))
profound=profoundProFound(image, magzero=30, verbose=TRUE, plot=TRUE)
#You can check to see if the final objects mask is aggressive enough. Notice the halos
#surrounding bright sources when just using the objects mask.
temp=image$imDat
temp[profound$objects>0]=0
magimage(temp)
temp=image$imDat
temp[profound$objects_redo>0]=0
magimage(temp)
magplot(profound$segstats[,c("R50","SB_N90")], log='x', grid=TRUE)
magplot(profound$segstats[,c("R50","SB_N90")], log='x', grid=TRUE)
magplot(profound$segstats[,c("flux","origfrac")], log='x', grid=TRUE)
# }
Run the code above in your browser using DataLab