Basic Testing of the Duchamp Source Finder
aa r X i v : . [ a s t r o - ph . I M ] D ec Basic Testing of the
Duchamp
Source Finder
T. Westmeier A , C , A. Popping A , and P. Serra B A International Centre for Radio Astronomy Research, M468, The University of Western Australia,35 Stirling Highway, Crawley WA 6009, Australia B ASTRON, Postbus 2, 7990 AA Dwingeloo, the Netherlands C Email: [email protected]
Abstract:
This paper presents and discusses the results of basic source finding tests in three dimensions(using spectroscopic data cubes) with
Duchamp , the standard source finder for the Australian SKAPathfinder. For this purpose, we generated different sets of unresolved and extended H i model sources.These models were then fed into Duchamp , using a range of different parameters and methods providedby the software. The main aim of the tests was to study the performance of
Duchamp on sources withdifferent parameters and morphologies and assess the accuracy of
Duchamp ’s source parametrisation.Overall, we find
Duchamp to be a powerful source finder capable of reliably detecting sources down tolow signal-to-noise ratios and accurately measuring their position and velocity. In the presence of noisein the data,
Duchamp ’s measurements of basic source parameters, such as spectral line width andintegrated flux, are affected by systematic errors. These errors are a consequence of the effect of noiseon the specific algorithms used by
Duchamp for measuring source parameters in combination with thefact that the software only takes into account pixels above a given flux threshold and hence misses partof the flux. In scientific applications of
Duchamp these systematic errors would have to be correctedfor. Alternatively,
Duchamp could be used as a source finder only, and source parametrisation couldbe done in a second step using more sophisticated parametrisation algorithms.
Keywords: methods: data analysis
With the advent of the
Square Kilometre Array (SKA;Dewdney et al. 2009) and its precursors and pathfind-ers, including the
Australian SKA Pathfinder (ASKAP;DeBoer et al. 2009), the
Karoo Array Telescope (Meer-KAT; Jonas 2009), and the
Aperture Tile In Focus (APERTIF; Oosterloo et al. 2009), the prospect ofdeep radio continuum and H i surveys of large areason the sky demands for new strategies in the areas ofdata reduction and analysis, given the sheer volume ofthe expected data streams, in particular for spectro-scopic surveys.Of particular importance is the automatic and ac-curate identification and parametrisation of sourceswith high completeness and reliability. Due to thelarge data volumes to be searched, source finding algo-rithms must be fully automated, and the once commonpractice of source finding ‘by eye’ will no longer be fea-sible. Moreover, accurate source parametrisation algo-rithms need to be developed to generate reliable sourcecatalogues free of systematic errors, as otherwise theintegrity of scientific results based on the survey datacould be compromised.In this paper we will take a closer look at the Duchamp source finder (Whiting 2011a, 2012). Du-champ has been developed by Matthew Whiting atCSIRO as a general-purpose source finder for three-dimensional data cubes as well as two-dimensional im- Duchamp ages and will serve as the default source finder in theprocessing of data from the ASKAP survey scienceprojects. The software identifies sources by searchingfor regions of emission above a specified flux thresh-old. To improve its performance,
Duchamp offers sev-eral methods of preconditioning and filtering of theinput data, including spatial and spectral smoothingas well as reconstruction of the entire image or datacube with the help of wavelets. In addition to sourcefinding,
Duchamp provides the user with basic sourceparametrisation, including the measurement of posi-tion, size, radial velocity, line width, and integratedflux of a source. More information about the capa-bilities of the software is available from the
Duchamp
User Guide (Whiting 2011b). A brief overview of
Du-champ ’s basic functionality is provided in Section 3.So far,
Duchamp ’s source finding and parametrisa-tion capabilities have never been systematically testedon a large set of artificial sources with well-defined pa-rameters. The aim of this paper is to bridge this gapby thoroughly testing the performance of
Duchamp on sets of artificial point sources and galaxy models aswell as a data set containing real galaxies and telescopenoise. The tests were originally motivated by the needto identify suitable source finding algorithms for the
Widefield ASKAP L-band Legacy All-sky Blind Sur-vey (WALLABY; Koribalski & Staveley-Smith 2009), one of the large, extragalactic ASKAP survey science Publications of the Astronomical Society of Australia
Table 1: Summary of the parameters used to gen-erate the visibility data set and noise image for thepoint source models.Parameter (visibility) Value UnitNumber of antennas 36System temperature 50 KDeclination − ◦ Total integration time 8 hHour angle range ± .
42 GHzChannel width 18 .
31 kHz3 .
86 km s − Parameter (image) Value UnitFinal image size 31 ×
31 pxField diameter 5 arcminPixel size 10 arcsecRobustness 0Gaussian uv taper 7 .
28 k λ .
54 km rms noise 1 .
95 mJySynthesised beammajor axis 27 . . . ◦ projects currently in preparation (Westmeier & John-ston 2010). Hence, the tests presented here will focuson the detection of compact and extended H i sources,in particular galaxies, in three-dimensional data cubeswith ASKAP characteristics.However, we believe that the results and conclu-sions presented in this paper will be of interest notonly to those involved in SKA precursor science, butto a larger community of astronomers interested in theautomatic detection and parametrisation of sources intheir data sets, regardless of the wavelength range in-volved. For a comparison of Duchamp ’s performancewith that of other source finding algorithms we referto the paper by Popping et al. (2011) in this issue.This paper is organised as follows: In Section 2 wesummarise the source finding strategies of other largeH i surveys in the past, followed by a brief overview ofthe Duchamp source finder in Section 3. In Section 4we present the outcome of our test of
Duchamp onpoint source models with simple Gaussian line pro-files. Section 5 describes our testing of
Duchamp onmodels of disc galaxies with varying physical parame-ters. In Section 6 we apply
Duchamp to a data cubecontaining real galaxies and genuine noise extractedfrom radio observations. A discussion of our results ispresented in Section 7. Finally, Section 8 summarisesour main results and conclusions.
Some of the previous, large H i surveys, including theH i Parkes All-Sky Survey (HIPASS; Barnes et al. 2001),the H i Jodrell All-Sky Survey (HIJASS; Lang et al.2003), and the Arecibo Legacy Fast ALFA survey (AL-FALFA; Giovanelli et al. 2005), already had to dealwith the issue of (semi-)automatic source detection.In the case of HIPASS, two different source find-ers were used and the results combined to maximisecompleteness (Meyer et al. 2004). The first algorithm, multifind (Kilborn 2001), used a simple 4 σ flux thresh-old combined with smoothing of the data cube on dif-ferent scales. The second algorithm, tophat , detectedsources in the spectral domain by convolving each spec-trum in the data cube with a top-hat function of vary-ing width. Neither of the two algorithms alone man-aged to detect more than 90% of the final, combinedsource list. The two algorithms combined producedabout 140,000 unique detections, all of which were in-spected by eye to remove potential false detections.The final HIPASS catalogue included 4315 sources (Meyeret al. 2004), resulting in an overall reliability of the au-tomatic source finding algorithms of only about 3%.In the case of HIJASS, again two different meth-ods were used and the results combined to improvecompleteness (Lang et al. 2003). The first methodsimply involved searching the cubes by eye to extractpotential sources. The second algorithm, polyfind ,first searched for signals above a given threshold ina smoothed version of the data cube and then ran aseries of matched filters over the detected signals todecide whether a signal was likely genuine or false. Asin the case of HIPASS, the source list produced by theautomatic source finding routine was inspected by eyeto further reject potential false detections. The posi-tions of uncertain detections were re-observed to eitherconfirm of refute them.For the ALFALFA survey, a matched-filtering tech-nique was applied to the data in the spectral domain(Saintonge 2007). The data were convolved with a setof template functions created by combining the firsttwo symmetrical Hermite functions, Ψ ( x ) and Ψ ( x ).The resulting templates range from simple Gaussianprofiles for narrow signals to double-peaked profiles forbroader signals, covering the range of spectral shapesexpected from H i observations of real galaxies. In testson 1500 simulated galaxies, 100% reliability and about70% completeness are achieved at an integrated signal-to-noise ratio of S/N ≈
6, while the 90% completenesslevel is exceeded at
S/N & Duchamp
Source Finder
Duchamp has been implemented as a general-purposesource finder for three-dimensional spectral-line datacubes with two spatial axes and one frequency (or In her calculation of
S/N , Saintonge (2007) makes theimplicit assumption that the sources are spatially unre-solved. ww.publish.csiro.au/journals/pasa velocity) axis, although the software can also oper-ate on one- and two-dimensional data sets (Whiting2011b). Duchamp finds sources by applying a sim-ple flux threshold to the data cube, specified by the threshold or snrCut keyword, and searching for pix-els above that threshold. In a second step, the soft-ware attempts to merge detections into sources underspecific circumstances that can be controlled by theuser. One option is to simply merge adjacent pix-els ( flagAdjacent keyword). Alternatively, a max-imum spatial and spectral separation can be speci-fied for the merging of detected pixels into sources( threshSpatial and threshVelocity keywords, respec-tively). Once detected, sources can be “grown” toa flux level below the actual detection threshold, us-ing the flagGrowth and growthThreshold / growthCut keywords.Basic removal of false detections is achieved by re-quiring that sources comprise a minimum number ofcontiguous spatial pixels and spectral channels, us-ing the minPix and minChannels keywords, respec-tively. To improve the reliability of the source find-ing even further, Duchamp offers a powerful methodof reconstructing the entire input data cube with thehelp of wavelets. Source finding is then performedon the reconstructed cube instead of the original in-put cube. Reconstruction can either be carried outin all three dimensions of the cube, or in the spa-tial (two-dimensional) or spectral (one-dimensional)domain only.
Duchamp uses the so-called ‘`a trous’ wavelet re-construction method (Starck & Murtagh 2002). First,the input data set is convolved with a specific waveletfilter function (three different functions are offered tothe user by
Duchamp ). The difference between theconvolved data set and the original data set is thenadded to the reconstructed cube. Next, the scale ofthe filter function is doubled and the procedure re-peated, using the convolved array as the new inputdata set. Once the user-specified maximum filter scaleis reached, the final convolved data set is added to thereconstructed cube, and source finding on the recon-structed data set commences.The ‘`a trous’ wavelet reconstruction of the datacube offers a powerful method of enhancing
Duchamp ’ssource finding capabilities. First of all, the user canselect the minimum ( scaleMin keyword) and maxi-mum ( scaleMax keyword) filter scales to be used inthe reconstruction, providing efficient suppression ofsmall-scale and large-scale artefacts in the data, suchas noise peaks, baseline ripples, or radio-frequency in-terference. Furthermore, the user can specify an ad-ditional threshold ( snrRecon keyword) to be appliedwhen adding wavelet components to the reconstructeddata cube, thereby reducing even further the numberof spurious signals in the data cube.In comparison to simple data thresholding, the ‘`atrous’ wavelet reconstruction method will greatly in-crease the completeness and reliability of
Duchamp ’ssource finding procedure, and hence the method hasbeen applied in all source finding tests presented inthis paper.
Final data cubesAdd random portion ofnoise cube to each cubewith beam modelConvolve each cubepoint sources with imgenCreate 1024 cubes ofGenerate visibiliydata with uvgenvisibilities with invertFourier−transformnoise cube beam model
Figure 1: Outline of the procedure used to cre-ate the model data cubes of point sources with
Miriad . For our first test of
Duchamp we generated modelsof 1024 point sources with simple Gaussian spectralline profiles. This will allow us to assess the funda-mental performance of
Duchamp under ideal condi-tions and to investigate the accuracy of the software’ssource parametrisation algorithms. Point sources withGaussian profiles are ideal for this test because—as aconsequence of their simple morphology—their physi-cal parameters can be exactly defined and calculatedto serve as a benchmark for
Duchamp ’s parametrisa-tion.In order to create the model data set, the
Miriad (Sault, Teuben, & Wright 1995) task uvgen was em-ployed to generate visibility data of Gaussian noise at afrequency of 1 . Miriad ’s task invert to generate a noise image of600 ×
600 pixels and 31 spectral channels with char-acteristics similar to WALLABY (again, see Table 1for details). The rms noise level in this image is σ =1 .
95 mJy which is only slightly higher than the 1 . Miriad task imgen was used to create 1024 data cubeseach of which has a size of 31 ×
31 pixels and 31 spectralchannels and contains a single point source with Gaus-sian spectral line profile in the centre. Each source wasrandomly assigned a peak flux in the range of 1 to 20 σ ,resulting in an average of about 54 sources per 1 σ in- Publications of the Astronomical Society of Australia
HPBW
Figure 2: Example of a point source model generated for testing
Duchamp . The left-hand panel showsa single-channel map of the data cube (at the systemic velocity of the source), and the right-hand paneldepicts the spectrum at the source position. The circle in the map illustrates the half-power beam width. terval. Spectral line widths (FWHM) range from 0 . . . − , resulting in a density of about27 sources per 1 km s − line width interval. Whilein reality sources with H i line widths of as small as0 . − will not exist, the reason for including suchnarrow lines in our test is to study the performance of Duchamp on sources that are spectrally unresolved,irrespective of the absolute line width.Each of the 1024 cubes was convolved with thebeam model produced by invert . Next, a random por-tion of 31 ×
31 pixels of the original noise cube was se-lected and added to each convolved image to create thefinal images used for testing
Duchamp . To facilitatecorrect integrated flux measurements, we added infor-mation on the synthesised beam to the image header.The entire procedure is outlined in Figure 1. An ex-ample image and spectrum of one of the point sourcemodels is shown in Figure 2.
Duchamp
Next, we ran
Duchamp (version 1.1.8) on the datacubes. In order to find out which combination of con-trol parameters provided the best performance in termsof completeness and reliability, we first ran
Duchamp several times with different flux thresholds and min-imum wavelet scales to test the performance of eachset of parameters. An overview of the different com-pleteness and reliability levels achieved in these runsas a function of integrated flux of the source is shownin Figure 3.We then selected the best set of control parametersfor the analysis presented in this section. In this best-performing run (number 5 in Figure 3 and Table 2) weused a 1 . σ flux threshold equivalent to 2 . Duchamp ’s ‘`a trous’ waveletreconstruction. We employed a full three-dimensionalwavelet reconstruction with a minimum wavelet scaleof 2 (i.e. the smallest scales were excluded to suppressnoise in the reconstructed cube) and a flux threshold of 3 σ for wavelet components to be included in the re-constructed cube. In addition, we required sources tocover a minimum of 5 contiguous pixels in the imagedomain and 3 contiguous spectral channels above thedetection threshold to be included in the final sourcecatalogue. This will further reduce the number of spu-rious detections. The Duchamp input parameters ex-plicitly set in the parameter file are listed in Table 3.The 1024 output parameter files generated by
Du-champ were concatenated, and those source entrieswhose positions were within ± Duchamp to processa certain amount of data. Firstly, the performance of
Duchamp strongly depends on the exact choice of in-put parameters, including detection threshold, waveletreconstruction choices, or settings related to mergingand discarding of initial detections. Three-dimensionalwavelet reconstruction of the input data cube, for ex-ample, is particularly computationally expensive. Sec-ondly, the running time of
Duchamp on a particulardata cube will depend on a large number of details,including the number of sources in the cube, their sizeand morphology, and in principle even the number den-sity of sources in the cube. Thirdly, recent updates ofthe software have resulted in a significant improvementof
Duchamp ’s performance, in particular compared toversion 1.1.8 used for part of the testing presented inthis paper.To get a basic idea of the impact of the aforemen-tioned parameters on the running time of
Duchamp we performed a few simple tests on a standard laptopcomputer with a state-of-the-art, dual-core 2 . Duchamp several times with different parameters on our arti-ficial noise data cube of 600 ×
600 spatial pixels and ww.publish.csiro.au/journals/pasa C o m p l e t e n e ss Completeness (Integrated Flux)Duchamp 1 Rel=63.0%Duchamp 2 Rel=97.0%Duchamp 3 Rel=99.1%Duchamp 4 Rel=99.2%Duchamp 5 Rel=63.2%Duchamp 6 Rel=97.0%Duchamp 7 Rel=98.6%Duchamp 8 Rel=6.3%
Figure 3: Completeness as a function of integratedflux for different tests of
Duchamp with varyingcontrol parameters. The parameters employed inthe different runs are listed in Table 2. The overallreliability for each run is listed in the legend.
31 spectral channels without any sources in it. Usinga 5 σ detection threshold, Duchamp takes about 0 .
64 sof CPU time to run, producing no detections. Whenperforming a one-dimensional wavelet reconstructionin the spectral dimension prior to source finding, therunning time increases by a factor of 30 to about 19 s.Full three-dimensional wavelet reconstruction is slowerby a factor of 120, requiring about 77 s to complete.As mentioned before, these numbers strongly de-pend on the number and nature of sources present inthe cube. Decreasing the flux threshold to 3 σ withoutwavelet reconstruction results in 966 detections (allof which are noise peaks) and increases the runningtime of Duchamp to about 3 . . . Table 2: Relevant input parameters for the differ-ent test runs of
Duchamp in order to find the opti-mal set of control parameters (see Figure 3). Thebest-performing parameter set (run 5) was thenused for the analysis presented in this paper.Run threshold scaleMin1 1 . . . . . . . . dicating that an increase in cube size does not translateinto a proportional increase in processing time whendealing with wavelet reconstruction.In summary, the time Duchamp needs to processa data cube is a complicated function of not only themachine specifications (e.g. CPU, memory, data trans-fer speed), but also the input parameters (e.g. fluxthreshold, wavelet reconstruction) and properties ofthe data set concerned (e.g. cube dimensions, numberof sources). Hence, running times are almost impos-sible to predict and may have to be determined ex-perimentally on a case-by-case basis. Instead of ask-ing whether
Duchamp is fast enough for a particularproblem, the user would have to determine the opti-mal set of conditions that would allow processing ofthe data in a given period of time. An alternative op-tion would be to separate the problem into multiple,parallel processes.
Two of the most important parameters in the charac-terisation of source finder performance are complete-ness and reliability. Completeness is defined as thenumber of genuine detections divided by the true num-ber of sources present in the data. Completeness caneither be calculated for the entire sample or more sen-sibly for a subset, e.g. for sources within a certain pa-rameter range. Reliability is defined as the number ofgenuine detections divided by the total number of de-tections produced by the source finder. Reliability canonly be calculated for the entire sample of sources andnot for a subset of sources within a certain parameterrange, because false detections do not possess physicalparameters as such. Alternatively, the parameters de-rived by the parametrisation algorithm of the sourcefinder can be used to derive reliability as a functionof different source parameters, but it is important tonote that for genuine sources those parameters can beaffected by systematic errors and do not necessarilycorrespond to the original source parameters.The ideal source finder would produce a complete-ness and reliability of 100%. In reality, however, we
Table 3:
Duchamp input parameters (Whiting2011b) explicitly set in the input parameter filefor point source models. The default values of
Duchamp were used for all other parameters.Parameter Value Commentthreshold 0.0029265 1.5 × rms minPix 5minChannels 3flagAdjacent trueflagATrous true Wavelet reconstr.reconDim 3 in 3 dimensionssnrRecon 3scaleMin 2 Publications of the Astronomical Society of Australia C o m p l e t e n e ss ( % ) C o m p l e t e n e ss ( % ) C o m p l e t e n e ss ( % ) Figure 4:
Top panel:
Completeness of the pointsource models as a function of true peak signal-to-noise ratio in bins of 1 σ . Middle panel:
Same,but as a function of true integrated flux in binsof 0 . − . Bottom panel:
Same, but as afunction of true line width (FWHM) in bins of2 . − . will have to find a compromise between good complete-ness and good reliability. In the case of Duchamp , forexample, decreasing the flux threshold for detectionswill lead to an increase in completeness, but at the costof lower reliability.In our test of
Duchamp on the set of 1024 pointsource models the software finds a total of 1103 sourcesof which 850 are genuine detections. The remaining253 detections are false positives due to strong noise R e li a b ilit y ( % ) R e li a b ilit y ( % ) R e li a b ilit y ( % ) Figure 5:
Top panel:
Reliability of the point sourcemodels as a function of measured peak signal-to-noise ratio in bins of 1 σ . Middle panel:
Same, butas a function of measured integrated flux in binsof 0 .
05 Jy km s − . Bottom panel:
Same, but as afunction of measured line width ( w ) in bins of2 . − . peaks in the data cube. These numbers translate intoan overall completeness of 83.0% and an overall relia-bility of 77.1%. Note that these numbers differ slightly from the onesquoted for run 5 in Figure 3 because a different realisationof the model was used in the initial tests. Reliability val-ues will generally depend on the characteristics of the datacube under consideration (e.g. the size of the cube) and aretherefore difficult to assess and compare. ww.publish.csiro.au/journals/pasa -10-50510-10-50510 Right ascension error (arcsec) D ec li n a ti on e rr o r ( a r c s ec ) P o s iti on e rr o r ( a r c s ec ) Figure 6:
Left-hand panel:
Position error of the point source models in right ascension and declination.
Right-hand panel:
Mean position error (black data points) and corresponding standard deviation (errorbars) as a function of true peak signal-to-noise ratio in 1 σ bins. Completeness as a function of peak signal-to-noiseratio is plotted in the top panel of Figure 4. The detec-tion list produced by
Duchamp is complete down toa peak flux level of F peak ≈ σ , but below that levelcompleteness decreases to below 50% at about 3 σ . Thecompleteness curve shows a much steeper rise whenplotted against integrated flux instead of peak flux(middle panel of Figure 4). The 100% completenesslevel is reached at F int ≈ . − , correspond-ing to an H i mass of about 7 × M ⊙ at a dis-tance of 1 Mpc, or 7 × M ⊙ at 100 Mpc, for theexpected 8-hour integration per pointing of the WAL-LABY project on ASKAP. Below that flux level thereis a sharp drop in completeness.The bottom panel of Figure 4 shows completenessplotted as a function of true line width (FWHM), ir-respective of the peak flux and integrated flux of asource. Over most of the covered line width rangethe completeness remains constant at approximately90%, but it gradually decreases to about 50% belowline widths of 10 km s − . This decrease is presumablythe result of the ‘`a trous’ wavelet reconstruction of thedata cube. By ignoring the smallest wavelet scales inthe reconstruction we suppress the detection of noisepeaks, but at the same time we are also less sensitiveto genuine sources with narrow spectral lines.Reliability as a function of measured peak signal-to-noise ratio, measured integrated flux, and measuredline width ( w ) is plotted in the top, middle, and bot-tom panels of Figure 5. Duchamp achieves 100% reli-ability at a peak signal-to-noise ration of about 5 andan integrated flux level of approximately 0 . − .Reliabilities range between about 80% to 100% overmost of the covered line width range, but drop signifi-cantly for sources with narrow lines of w .
15 km s − due to the increasing number of false detections asso-ciated with noise peaks. For line widths of less thanabout 5 km s − the reliability increases sharply, be-cause Duchamp effectively filters narrow signals causedby noise peaks through wavelet filtering and minimumchannel requirements. However, as discussed earlier, reliability calcula-tions are very difficult to assess and should be ap-proached with great caution. First of all, reliability canonly be specified as a function of measured source pa-rameters because false detections do not have genuinephysical parameters. Any errors in a source finder’sparametrisation will therefore affect the calculated re-liability curves. Secondly, the actual reliability num-bers are entirely meaningless in the case of model dataas they depend on how the sources were distributedacross the model data cube. Increasing the volumeof the cube (without increasing the number of sourcestherein) will result in lower reliabilities as the fractionof false detections increases. Consequently, reliabilitiescan only be compared on a relative scale, e.g. whentesting different source finders on the same data set todetermine which algorithm performs best.
The resulting position errors are plotted in the left-hand panel of Figure 6.
Duchamp does an excellentjob in determining accurate source positions, with amean position error of 0 . ± . . ± . σ , is shown inthe right-hand panel of Figure 6. For bright sources of F peak ≈ σ the mean position error is approximately1 arcsec, increasing to about 5 arcsec for F peak ≈ σ .These numbers correspond to only about 4% and 19%,respectively, of the FWHM of the synthesised beam.Two limitations should be noted at this point. Firstof all, in our models the source was always placed ex-actly on the central pixel of the data cube. We didnot explicitly test placement of sources at positionsin between the grid points of the cube, which—in thecase of point sources—could result in reduced detec-tion rates and less accurate source positions. Secondly,as with other source parameters, source positions will Publications of the Astronomical Society of Australia -4 -2 0 2 401020304050 Velocity error (km/s) N u m b e r o f s ou r ce s V e l o c it y e rr o r ( k m / s ) Figure 7:
Left-hand panel:
Histogram of radial velocity errors (black curve) of the point source modelsin bins of 0 . − . The red, dashed curve is the result of a Gaussian fit to the histogram. Right-handpanel:
Standard deviation of the velocity error of the sources as a function of true peak signal-to-noiseratio in bins of 1 σ . be inaccurate in cases where two or more sources areconfused. An overall histogram of radial velocity errors derivedfrom the
Duchamp run is shown in the left-hand panelof Figure 7. As expected, velocity errors have an ap-proximately Gaussian distribution centred on zero. Themean velocity error of all sources is 0 . ± . − .The red, dashed curve in Figure 7 shows the result ofa Gaussian fit to the histogram. While the overall dis-tribution of velocity errors follows the fitted Gaussianfunction, there are a few significant deviations, namelya somewhat higher and sharper peak in the centre(which is slightly shifted into the negative range) andconspicuous ‘wings’ between 2 and 3 km s − (both pos-itive and negative) where source counts are systemat-ically too high with respect to the fit. The FWHMof the fitted Gaussian is 1 . ± .
04 km s − , and thecentroid is − . ± .
017 km s − which deviates fromzero by about 1 . σ , reflecting the aforementioned neg-ative offset of the peak of the histogram. These devi-ations from a pure Gaussian distribution are possiblycaused by digitisation effects related to the segmenta-tion of the frequency axis into discrete bins of 18 . .
86 km s − .The standard deviation of the radial velocity erroras a function of peak signal-to-noise ratio in 1 σ binsis shown in the right-hand panel of Figure 7. As ex-pected, the standard deviation from the mean (whichis essentially zero for all peak flux intervals) increaseswith decreasing peak flux. While for bright sources of F peak ≈ σ the standard deviation is below 1 km s − ,it increases to almost 6 km s − for faint sources nearthe 3 σ level. Figure 8 shows the ratio of measured line width ver-sus true line width (FWHM of the original Gaussian model) as a function of peak signal-to-noise ratio inbins of 1 σ . Duchamp determines three different typesof line width: w is the full width at 50% of the peakflux, w is the full width at 20% of the peak flux,and w vel is the full detected line width of the source,i.e. the width across all channels with detected flux.For a Gaussian line, w is equivalent to the FWHM,and the ratio of FWHM /w should therefore be 1.The relation between w and w in the case of aGaussian line is given by the constant factor of w w = 1 . . (1)Finally, the relation between w vel and w , again as-suming a Gaussian line profile, is defined via w vel w = (cid:20) log (cid:18) F thr F peak (cid:19)(cid:21) (2)where F thr = n × σ is the flux threshold used in the cal-culation of w vel . These theoretical relations are plot-ted in Figure 8 as the dashed lines for w (black), w (red), and w vel (blue; for F thr = 1 . σ ). Duchamp ’s measurement of w (black data points)is in excellent agreement with the expectation (black,dashed line) over a wide range of peak signal-to-noiseratios. Only for faint sources of F peak < σ are themeasured line widths on average slightly smaller thanthe true widths, but by no more than about 10 to 15%.In contrast, Duchamp ’s measurements of w and w vel (red and blue data points, respectively) are sys-tematically too large over most of the covered rangeof signal-to-noise ratio as compared to the theoreticalexpectation (red and blue dashed lines, respectively).Only for faint sources of F peak . σ do the measured w fall short of the theoretical ones. This result sug-gests that w is the most accurate measurement ofline width provided by Duchamp and should be usedinstead of w and w vel for the characterisation of as-tronomical sources. However, both w and w sys-tematically fall short of the true line width for faintsources below F peak ≈ σ . ww.publish.csiro.au/journals/pasa L i n e w i d t h r a ti o Figure 8: Ratio of measured versus true line widthfor the point source models as a function of truepeak signal-to-noise ratio in bins of 1 σ . The blackdata points show w , the red data points w , andthe blue data points w vel (for a 1 . σ flux thresh-old), all of which have been divided by the originalFWHM of the Gaussian line. The correspondingtheoretical expectations for a Gaussian line profileare shown as the black, red, and blue dashed lines. The ratio of recovered versus true peak flux of themodel point sources is plotted in the left-hand panel ofFigure 9 as a function of peak signal-to-noise ratio in1 σ bins. The dashed and dotted red lines indicate thetheoretical ± σ and ± σ envelopes, respectively. Theright-hand panel shows the same figure, but as a func-tion of integrated flux in bins of 0 . − . Forbright sources of F peak & σ Duchamp accuratelyrecovers the peak flux of the sources, although thereis the general tendency of measured peak fluxes be-ing slightly too high on average. For fainter sources of F peak . σ there is a strong deviation, with measuredfluxes being systematically too high by a significantfactor. This is generally due to faint sources beingmore likely to be detected when their maximum coin-cides with a positive noise peak, whereas faint sourcessitting on top of a negative noise peak will likely re-main undetected, thereby creating a strong bias in themeasurement of peak fluxes.Even for high peak signal-to-noise ratios the peakfluxes measured by Duchamp tend to be slightly toolarge.
Duchamp determines the peak flux of a sourceby simply selecting the data element with the highestflux encountered. As mentioned before, this method isbiased towards selecting data elements that have beenaffected by positive noise peaks. In sources with broadspectral signals there is a higher probability of findinga positive noise signal in one of the channels near thepeak of the line that increases the signal beyond the ac-tual line peak. This is a result of the source being wellresolved in the spectral domain. Hence, peak fluxesmeasured by
Duchamp will generally be too high ir-respective of source brightness as long as the sourceis spectrally (or spatially) resolved. For very bright sources, however, the relative error will be negligible.
The ratio of measured versus true integrated flux ofthe model point sources as a function of peak signal-to-noise ratio in bins of 1 σ is presented in the left-hand panel of Figure 10. The right-hand panel showsthe same figure, but as a function of integrated flux inbins of 0 . − . Apparently, Duchamp ’s mea-surement of the integrated flux of a source is system-atically too low by a significant factor. Even for thebrightest sources of F peak ≈ σ only about 90% ofthe true flux is recovered by Duchamp , and that figuredrops to well below 50% for faint sources of F peak < σ .This issue is likely caused by the fact that Duchamp only considers pixels above the detection threshold whencalculating the integrated flux. Pixels below the thresh-old, while potentially contributing significantly to theoverall flux of a source, are not included in the summa-tion carried out by
Duchamp , resulting in integratedfluxes being systematically too small.In order to study the expected decrease in the inte-grated flux measurement, let us assume a point sourcewith Gaussian line profile being observed with a tele-scope with radially symmetric Gaussian point spreadfunction (PSF), F ( x, y, v ) = F peak exp (cid:18) − x + y σ − v σ v (cid:19) , (3)with amplitude, F peak , velocity dispersion, σ v , andPSF size, σ PSF . The integrated flux measurement canthen be considered as the integral under the three-dimensional Gaussian brightness profile across the fre-quency/velocity range, ± v , and the spatial range, ± x and ± y , over which the flux of the line is above thedetection threshold, thus F int = x Z − x y Z − y v Z − v F ( x, y, v ) d x d y d v (4)= F peak (2 π ) / σ σ v erf (cid:18) x √ σ PSF (cid:19) × erf (cid:18) y √ σ PSF (cid:19) erf (cid:18) v √ σ v (cid:19) , (5)where erf( x ) is the error function. Inserting the appro-priate integration limits and then dividing Equation 5by the total flux (i.e. integrated over ±∞ ) leads to atheoretical integrated flux ratio of F int F tot = h erf (cid:16)p − ln ( F thr /F peak ) (cid:17)i (6)with a flux threshold of F thr = n × σ .The resulting theoretical integrated flux accordingto Equation 6, assuming a 1 . σ threshold, is shownas the solid red curve in Figure 10. The integratedfluxes measured by Duchamp are only slightly belowwhat one would expect from a simple integration overa three-dimensional Gaussian. A fit to the data pointsinstead yields an effective flux threshold of 2 . σ , shownas the dotted red curve in Figure 10, which is slightly Publications of the Astronomical Society of Australia P ea k f l ux r a ti o P ea k f l ux r a ti o Figure 9:
Left-hand panel:
Ratio of measured versus true peak flux (black data points) and correspondingstandard deviation (error bars) of the model point sources as a function of true peak signal-to-noise ratioin 1 σ bins. The dashed and dotted red lines indicate the theoretical ± σ and ± σ envelopes, respectively. Right-hand panel:
Same, but as a function of true integrated flux in bins of 0 . − . I n t e g r a t e d f l ux r a ti o I n t e g r a t e d f l ux r a ti o Figure 10:
Left-hand panel:
Ratio of measured versus true integrated flux (black data points) and corre-sponding standard deviation (error bars) for the model point sources as a function of true peak signal-to-noise ratio in bins of 1 σ . The solid red curve shows the theoretical expectation for the 1 . σ flux thresholdused in our test. The dotted red curve shows the best fit to the data points, corresponding to an effec-tive flux threshold of 2 . σ . Right-hand panel:
Same, but as a function of true integrated flux in bins of0 . − . larger than the 1 . σ used when running Duchamp . Itis not quite clear why
Duchamp performs worse thanexpected. The discrepancy could be due to the factthat the software sums over discrete pixels whereas weassumed continuous integration in our mathematicalmodel. This will likely result in small differences, par-ticularly in those cases where the number of elementsacross the Gaussian profile is small. In our case, as weare dealing with point sources, this is certainly true forthe spatial dimension.In summary, integrated flux measurements pro-vided by
Duchamp are systematically too small andwill need to be corrected substantially to compensatefor the systematic offset.
In order to test the performance of
Duchamp on morerealistic, extended sources, we generated 1024 artificialH i models of galaxies with a wide range of parameters,using a programme written in C for direct manipula-tion of FITS data cubes. All galaxies were modelledas infinitely thin discs with varying inclination (0 ◦ to89 ◦ ), position angle (0 ◦ to 180 ◦ ), and rotation veloc-ity (20 to 300 km s − ). For any one galaxy, inclinationand position angle were considered to be constant overthe radial extent of the disc, while the rotation velocityincreases linearly from 0 to v rot between the centre and0 . .
65 km s − (equivalent to 2 . ww.publish.csiro.au/journals/pasa Duchamp . The left-hand panel shows thezeroth moment of the model, the middle panel shows the position-velocity diagram along the dashed, redline, and the right-hand panel depicts the integrated spectrum of the model galaxy. width of a spectral channel). The radial surface bright-ness profile was assumed to be Gaussian, too, resultingin an elliptical Gaussian brightness distribution on thesky.Next, we again generated an artificial ASKAP visi-bility data set of pure Gaussian noise at a frequency of1 . Miriad task uvgen with parametersas listed in Table 4. The visibility data were Fourier-transformed to create an image of the point spreadfunction and a noise data cube with 1601 × ×
32 galaxies and added to the noise cube to createthe final data cube of model galaxies for the testing of
Duchamp .The moment-zero map, position-velocity map, andintegrated spectrum of one of the model galaxies isshown in Figure 11 for illustration. As with the pointsources, all galaxies were centred on a pixel, althoughfor extended sources we do not expect any significanteffect from shifting the source centre with respect tothe pixel centre. Again, all sources are isolated, andwe did not attempt to test
Duchamp in a situationwhere source crowding occurs.It is important to note at this point that the re-sulting model galaxies, while exhibiting some of thespatial and spectral characteristics of real spiral galax-ies, have been simplified to a great extent, resulting inlimitations that need to be kept in mind when inter-preting the results presented in this section. Firstly,the assumption of an infinitely thin disc will result inunrealistic edge-on galaxies, with integrated fluxes aswell as individual spectral line widths across the discbeing too small. Secondly, parameters such as peakflux, angular size, or rotation velocity were all variedindependently of each other, resulting in unrealisticcombinations of galaxy parameters in some cases. Thepurpose of the models is to cover a vast parameter range of extended sources irrespective of whether thatentire range is populated by real galaxies. Even if discgalaxies with a certain combination of parameters donot exist, other objects, such as irregular galaxies orhigh-velocity clouds, could still cover those regions ofparameter space, and their exploration will thereforebe meaningful.
Table 4: Summary of the parameters used to gen-erate the visibility data set and noise image for thegalaxy models.Parameter (visibility) Value UnitNumber of antennas 36System temperature 50 KDeclination − ◦ Total integration time 8 hHour angle range ± .
42 GHzChannel width 18 .
31 kHz3 .
86 km s − Parameter (image) Value UnitFinal image size 1601 × . ◦ uv taper 7 .
28 k λ .
54 km rms noise 1 .
86 mJySynthesised beammajor axis 30 . . . ◦ Publications of the Astronomical Society of Australia
Table 5:
Duchamp input parameters explicitly set in the input parameter file for the galaxy models. Thedefault values of
Duchamp were used for all other parameters.Parameter Run 1 Run 2 Run 3 Commentthreshold 0.00186 0.00186 0.00186 1.0 × rms minPix 10 10 10minChannels 5 5 3flagAdjacent true true trueflagGrowth false true truegrowthThreshold – 0.00093 0.00093 0.5 × rms flagRejectBeforeMerge false true trueflagATrous true true true Wavelet reconstructionreconDim 3 3 3 in 3 dimensionssnrRecon 2 2 2scaleMin 3 3 3 Duchamp
We ran
Duchamp (version 1.1.12) on the model galaxycube several times with slightly different input param-eters to compare the performance. The different inputparameters explicitly set in the parameter file are listedand compared in Table 5. In all cases we employed a1 σ flux threshold, equivalent to about 1 . σ for wavelet components to be included in thereconstructed cube. The slightly larger minimum scaleas compared to the point source models is motivatedby the fact that we are now dealing with spatially andspectrally much more extended sources. In addition,we varied the number of contiguous spectral channelsrequired for detections and used Duchamp ’s growthcriterion in a few runs with a growth flux threshold of0 . σ . The latter method will grow detections to fluxlevels below the original detection threshold, resultingin more accurate source parametrisation. As it turnedout, the change from 5 to 3 consecutive spectral chan-nels for detections (run 2 versus 3) did not have anymajor impact on the results. Hence, only the resultsof runs 1 and 3 will be presented and discussed here.In order to compare the outcome of Duchamp withthe original input catalogue, we wrote a short Pythonscript that reads in and processes the different cata-logues. The script first reads in the
Duchamp outputcatalogue, the original model catalogue, and a spe-cial mask data cube marking pixels with emission inthe original model by assigning them a unique numbercharacteristic to each input source. The script thencycles through all the detections made by
Duchamp and decides for each detection whether it is genuine ornot by checking the value of the mask data cube at thesame position. If the detection is found to be genuine,the script will cycle through the original model cat-alogue to extract the actual input parameters of therespective source for comparison with the parametri-sation results of
Duchamp .At the end of this process we get a match of de-tected sources with original input sources, allowing usto calculate parameters such as completeness, reliabil- ity, and the fraction of sources being broken up intomultiple detections by
Duchamp . In addition, we areable to compare the original parameters of each sourcewith those determined by
Duchamp to test the per-formance of
Duchamp ’s parametrisation algorithms.
For run 1 (without growing of detections to flux levelsbelow the threshold), 436 out of 1063 detected sourcesare genuine, resulting in an overall reliability of 41%.As many original sources got broken up into multipledetections, only 194 of the 1024 input galaxies were de-tected, yielding an overall completeness of only 19%.There is a significant improvement for run 3 (withgrowing of detections to a flux level of 0 . σ ), where542 out of 1051 detected sources are genuine (relia-bility of 52%), but this time 521 of the 1024 inputgalaxies were detected, resulting in a much improvedoverall completeness of 51%.Completeness as a function of different galaxy pa-rameters is shown in Figure 12 for runs 1 and 3 (blackand red data points, respectively). As mentioned be-fore, run 1 resulted in very low completeness values oftypically only about 20% and no strong variation witheither the integrated flux of a source or its inclinationand rotation velocity. By growing detections to a fluxlevel of 0 . σ (run 3) we achieved much higher complete-ness levels over a large parameter range. 100% com-pleteness is achieved for sources of F int &
20 Jy km s − ,and completeness levels reach 50% at F int ≈ . − .The latter corresponds to an H i mass sensitivity of6 × M ⊙ at a distance of 1 Mpc, or 6 × M ⊙ at 100 Mpc, for the expected 8-hour integration perpointing of the WALLABY project on ASKAP.As shown in the middle and bottom panels of Fig-ure 12, there is a strong variation of completeness withboth inclination and rotation velocity of the galaxies.While face-on galaxies are on average detected at com-pleteness levels near 80%, Duchamp struggles to findedge-on galaxies, yielding average completeness levelsof only about 20% for galaxies with inclination angles ww.publish.csiro.au/journals/pasa C o m p l e t e n e ss ( % ) C o m p l e t e n e ss ( % ) C o m p l e t e n e ss ( % ) Figure 12: Completeness for the galaxy modelsas a function of true integrated flux in bins of2 . − (top panel), galaxy inclination inbins of 5 ◦ (middle panel), and rotation velocityin bins of 19 . − (bottom panel) for runs 1(black) and 3 (red). greater than 80 ◦ . This effect is caused by the com-bination of two separate effects. Firstly, as a resultof the limitations from our assumption of an infinitelythin disc, edge-on galaxies have typically lower inte-grated fluxes than face-on galaxies. Secondly, edge-on galaxies typically have a broader spectral signatureas a result of their higher projected rotation velocity,making it more difficult for Duchamp to pick up theirextended signal.The latter effect can also be seen in the bottom B r ea k - up o f s ou r ce s ( % ) B r ea k - up o f s ou r ce s ( % ) B r ea k - up o f s ou r ce s ( % ) Figure 13: Fraction of model galaxies being bro-ken up into two or more separate detections by
Duchamp as a function of true integrated flux inbins of 2 . − (top panel), galaxy inclina-tion in bins of 5 ◦ (middle panel), and rotation ve-locity in bins of 19 . − (bottom panel) forruns 1 (black) and 3 (red). panel of Figure 12, where completeness levels system-atically decrease as a function of increasing rotationvelocity of a galaxy, irrespective of its inclination orintegrated flux, confirming that on average Duchamp is more likely to pick up face-on galaxies with narrowspectral lines. It is important to note, however, that ata given distance galaxies with higher rotation velocitywill typically have a larger H i mass and are thereforemore likely to be detected than galaxies with lower Publications of the Astronomical Society of Australia P o s iti on e rr o r ( a r c s ec ) -40-2002040-40-2002040 Right ascension error (arcsec) D ec li n a ti on e rr o r ( a r c s ec ) Figure 14:
Left-hand panel:
Position errors (from run 3) for the model galaxies in right ascension anddeclination.
Right-hand panel:
Mean absolute position error (data points) and standard deviation (errorbars) as a function of true integrated flux in bins of 2 . − . rotation velocity at the same distance. Due to their rotation velocity, spiral galaxies often ex-hibit a large radial velocity gradient across their pro-jected disc on the sky, resulting in the characteristicdouble-horn profile of their integrated spectrum. This,however, can result in the two halves of a galaxy beingdetected as two separate sources by
Duchamp , in par-ticular in the case of faint, edge-on galaxies with largerotation velocities.In Figure 13 we have plotted the fraction of de-tected model galaxies that were broken up into two ormore separate detections by
Duchamp as a functionof true integrated flux (top panel), inclination (mid-dle panel), and rotation velocity (bottom panel). Forrun 1 (black data points) there is a very high frac-tion of multiple detections, typically about 60–80%,with no strong variation with either integrated flux ofthe galaxy or inclination and rotation velocity of thedisc. In total, 136 out of the 194 detected galaxies, or70 . Duchamp .Growing detections to the 0 . σ level (run 3, reddata points) results in a major improvement, with thenumber of multiple detections (in total 21 out of 521 de-tected galaxies, or 4 . F int . − does the fraction of multiple de-tections gradually increase up to about 10% at the lowend of the flux spectrum. Figure 13 also clearly showsthe expected increase in multiple detections for galax-ies of higher inclination ( i & ◦ ) and rotation velocity( v rot &
150 km s − ), which is the result of the double-horn profile becoming wider and more pronounced asthe radial velocity gradient in the plane of the sky in-creases.A similar case, although more difficult to assess, isthe detection of only one half of a galaxy (one horn of the double-horn profile), whereas the other half re-mains undetected. As there is only a single detection ofeach affected galaxy, such partial detections are muchmore difficult to identify. They should, however, resultin a significant offset of both the measured position andradial velocity of the detected source with respect tothe location of the originating galaxy.In the case of run 3, 62 out of 500 single detectionsshow velocity errors of more than 20 km s − , with 28even exceeding 150 km s − . The former corresponds toa fraction of 12 .
4% of all single detections. Similarly,62 out of 500 singly detected sources have a positionerror of more than 20 arcsec, which again correspondsto a fraction of 12 . These results suggest that, even when growing de-tections down to the 0 . σ level, there is a significantnumber of partial (approximately 66 sources) or multi-ple (21 sources) detections, corresponding to an overallfraction of about 16 .
7% of all genuine detections. Suchcases need to be identified in the output catalogue pro-duced by
Duchamp , as otherwise they will introducea significant bias in the measurement of source param-eters such as line width and H i mass. Identification ofbroken-up sources will be a very difficult task in prac-tice, as it may be impossible to decide whether twodetections are part of the same source or two sepa-rate sources in close proximity. While the growing ofdetections to lower flux levels can in principle reducethe fraction of sources being broken up, an undesirableside effect will be the potential merging of neighbour-ing sources, e.g. close galaxy pairs in group or clusterenvironments. The left-hand panel of Figure 14 shows a scatter plotof position errors for the model galaxies (based on There is no exact match between the 62 sources withlarge position error and the 62 sources with large velocityerror. A total of 66 sources fulfil either of the two criteria. ww.publish.csiro.au/journals/pasa -15 -10 -5 0 5 10 15020406080 Velocity error (km/s) N u m b e r Figure 15: Histogram of radial velocity errors(from run 3) for the model galaxies in bins of0 . − . run 3) in right ascension and declination. The meanposition errors in right ascension and declination are1 . ± . . ± . . ± . . ± . F int &
10 Jy km s − source positions are very accurate with typical errorsof about 2 . The mean velocity error (based on run 3) for the galaxymodels is − . ± . − . As in the case of sourceposition, the large standard deviation about the meanis caused by galaxies that are only partially detected.By including only sources with position errors of lessthan 15 arcsec in both right ascension and declinationand velocity errors of less than 20 km s − we can ex-clude such partial detections, resulting in a correctedmean radial velocity error of − . ± . − .A histogram of radial velocity errors for the galaxymodels is shown in Figure 15. As in the case of pointsources, the distribution is not exactly Gaussian. In-stead, there is a sharp peak near zero and an under-lying broad distribution of errors, in particular in thenegative range. Some of these non-Gaussian structures could again be the result of digitisation effects in con-junction with the spectral channel width of 3 .
86 km s − ,while we have no conclusive explanation for the notice-able asymmetry of the distribution. In order to estimate the original line width of the inputmodels, we calculated a ‘pseudo line width’ which bal-ances the intrinsic width of an individual line profilewith the overall, integrated line width resulting fromthe rotation velocity of the galaxy, thus w mod = q [2 v rot sin( i )] + w (7)where v rot is the rotation velocity of the model galaxy, i is the inclination of the disc, and w int = 22 . − is the intrinsic FWHM of the Gaussian spectral line ateach position across the galaxy.The left-hand panel of Figure 16 shows the meanratio of the measured line width, w , over the cal-culated ‘pseudo line width’, w mod , as a function oftrue integrated flux in bins of 2 . − (based onrun 3). Duchamp measures accurate line widths closeto the true value over a wide range of fluxes. The smalldeviation from the value of 1 can be easily explainedby the fact that w mod is just an approximation to theFWHM of the line profile. Only for fainter sources of F int . − does the line width ratio decreaseand the standard deviation increase significantly, in-dicating larger errors in Duchamp ’s measurement ofline width.In the right-hand panel of Figure 16 we have plot-ted the ratio of w /w mod as a function of w mod inbins of 50 km s − . While line width measurements forsources with narrow lines of w mod .
250 km s − are onaverage accurate, there is a systematic discrepancy forsources with broader lines, the line widths measured by Duchamp being systematically too small. The largestandard deviation suggests that this could have beencaused by cases in which only one half of the galaxywas detected, whereas the other half remained unde-tected, resulting in a significantly lower value of themeasured line width. Nevertheless, line width mea-surements for fully-detected sources should be accu-rate even if their line widths are large. This problemagain demonstrates the need to identify partially de-tected sources to avoid systematic errors that wouldaffect the scientific interpretation of the data.
The ratio of measured versus true integrated flux of themodel galaxies, based on run 3, is shown in Figure 17.Similar to our previous tests on point sources (see Fig-ure 10), the integrated flux measured by
Duchamp issystematically too low. For bright sources of F int ≈
20 Jy km s − a large fraction of approximately 95% ofthe flux is recovered, whereas this figure drops to be-low 60% for fainter sources of F int . − . Atthe same time, the scatter significantly increases, sug-gesting larger uncertainties (on a relative scale) in theflux measurement of faint sources. Publications of the Astronomical Society of Australia L i n e w i d t h r a ti o L i n e w i d t h r a ti o Figure 16: Ratio of measured ( w , from run 3) versus true line width for the model galaxies as a functionof true integrated flux in bins of 2 . − (left-hand panel) and true line width in bins of 50 km s − (right-hand panel). The error bars indicate the standard deviation about the mean. I n t e g r a t e d f l ux r a ti o Figure 17: Ratio of measured (from run 3) ver-sus true integrated flux of the model galaxies asa function of true integrated flux. The error barsindicate the standard deviation about the mean.The bin width decreases towards lower fluxes.
As discussed previously, the reason for the failureof
Duchamp to accurately determine the integratedflux of a source is that the software only sums over dataelements that are above the flux threshold and hencemisses some of the flux. Even the growth of detectionsdown to the 0 . σ level has not solved this fundamentalproblem, although the defect has become less severethan for the point source models without growing (seethe right-hand panel of Figure 10 for comparison). So far, we have tested
Duchamp on artificial sourcesembedded in perfectly Gaussian noise. While this isuseful to study the basic performance of the software,real observations will be more challenging for any source finder due to the more complex morphology of realsources and the presence of various artefacts in thedata, e.g. terrestrial and solar interference, spectralbaseline instabilities, or residual continuum emission.Unfortunately, it is not possible to simply test
Du-champ on a real H i data cube, because we would nothave a-priori knowledge of the sources in such a cubeand would not be able to assess which of the detec-tions made by Duchamp are genuine. A solution tothis problem would be to inject copies of real galax-ies into a real data cube of “pure” noise, i.e. a datacube extracted from telescopic observations that doesnot contain any H i sources above the noise level. Thismethod combines the advantages of artificial sourcemodels, where the source locations and parameters areexactly known, with those of real observations with re-alistic sources and artefacts.For this purpose, we generated a data cube con-taining real noise extracted from an observation withthe Westerbork Synthesis Radio Telescope (WSRT).We then added about 100 data cubes from the “West-erbork Observations of Neutral Hydrogen in Irregularand Spiral Galaxies” (WHISP) survey (Kamphuis, Si-jbring, & van Albada 1996; Swaters et al. 2002), eachcontaining one or more galaxies. The selected WHISPdata cube were artificially redshifted by scaling theirsize and flux level to match sources in a redshift rangeof 0 . . z . .
04, centred on the median redshift of0 .
03 expected for the WALLABY project (Koribalski& Staveley-Smith 2009). The procedure for creatingthe test data cube is explained in more detail by Serra,Jurek, & Fl¨oer (2011).The final test data cube has a size of 360 ×
360 spa-tial pixels and 1464 spectral channels. The pixel size of10 arcsec (with a synthesised beam width of 30 arcsec)and channel width of 18 . − ) were chosen to reflect the expected speci-fications of WALLABY. Figure 18 shows an exampleimage and spectra of two of the galaxies in the finalcube. As the locations and properties of the injectedgalaxies are well-known, we can directly compare them ww.publish.csiro.au/journals/pasa A B AB
Figure 18:
Left-hand panel:
Moment-zero map of a small region of the WSRT cube with injected WHISPgalaxies, showing two galaxies labelled A and B.
Right-hand panels:
Integrated spectra of the two galaxies. to the output of
Duchamp to assess performance in-dicators such as completeness and reliability.
Duchamp
We ran
Duchamp multiple times on the WSRT datacube with WHISP galaxies to probe different inputparameter settings of
Duchamp . A summary of theruns and parameters used is given in Table 6. Wemainly covered a wide range of flux thresholds between1 . . σ and tested one-dimensional (spectral do-main only) versus three-dimensional (spatial and spec-tral domain) wavelet reconstruction of the cube. Theoutput catalogue of each run was again cross-matchedwith the original source catalogue, using the Pythonscript described in Section 5.1.In order to obtain the original source catalogue, weran Duchamp once on the input model cube withoutnoise, using a very low detection threshold of well be-low the final noise level and no wavelet reconstruction.This resulted in a list of 100 sources against whichthe output catalogue provided by
Duchamp can bejudged. Since this method already introduces a strongbias in the catalogue of source parameters, we will onlyanalyse the completeness and reliability of
Duchamp ,but we shall not attempt to assess the parametrisationperformance of the software, because we do not havean exact source catalogue against which we would beable to assess the source parameters as measured by
Duchamp from the final test cube.
Completeness and reliability of the different runs of
Duchamp on the WSRT model cube with WHISPgalaxies are listed in Table 6 and displayed in Fig-ure 19 as a function of detection threshold. Generally,between about 40% to 60% of all galaxies in the cube were found by
Duchamp , while the overall reliabilityvaries strongly from about 10% to 100% depending ondetection threshold and wavelet reconstruction param-eters.We achieve better results for one-dimensional wave-let reconstruction (black and blue data points in Fig-ure 19) which generally yields higher completeness andreliability than three-dimensional wavelet reconstruc-tion (red data points). This is presumably due tothe small angular size of most galaxies in the modelcube; there is not much to gain from performing awavelet reconstruction in the spatial domain, whereasone-dimensional wavelet reconstruction in the frequencydomain yields much better results because most galax-ies are well-resolved and extended in frequency.In Figure 20 we plot completeness as a function ofintegrated flux for selected runs of
Duchamp . Above aflux of F int & − Duchamp consistently findsall sources irrespective of the input parameters chosen.At lower fluxes the different runs produce significantlydifferent results, with the one-dimensional wavelet re-construction (black and blue data points) generallyperforming better than the three-dimensional recon-struction (red data points), as noted before.The best-performing parameter set in terms of com-pleteness, run 7, produces a completeness of 50% at anintegrated flux of F int ≈ . − , correspondingto an H i mass of 1 . × M ⊙ at a distance of 1 Mpc,or 1 . × M ⊙ at 100 Mpc. This is worse than whatwe achieved for the point sources with Gaussian lineprofiles in Section 4, but significantly better than theoutcome for the model galaxies in Section 5. The rea-son for the better performance could be that the arti-ficially redshifted WHISP galaxies are generally muchmore compact than the model galaxies created for thetests in Section 5. As with any threshold-based sourcefinder, compact sources are easier to detect than ex-tended sources, even with prior wavelet reconstruction Publications of the Astronomical Society of Australia
Table 6: Relevant
Duchamp input parameters set in the different input parameter files for the modelbased on real WHISP galaxies and WSRT noise. The default values of
Duchamp were used for most of theother parameters. The last two rows list the overall completeness and reliability achieved by
Duchamp .Run 1 2 3 4 5 6 7 8 9 10 11threshold 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . or smoothing. As the spectral profiles of the WHISPgalaxies are generally broad and complex, performanceis worse than in the case of the point source models inSection 4 which had much simpler and narrower Gaus-sian lines.The overall reliability in the case of run 7 is verylow with only 10%. This figure, however, is the rawreliability achieved by Duchamp and can be substan-tially improved by filtering sources based on their mea-sured parameters. False detections are usually the re-sult of noise peaks being picked up by the source finder.A large fraction of these false noise detections will becharacterised by very low integrated fluxes and smallline widths, and often a simple cut in flux–line widthspace will remove more than 95% of false detectionswhile retaining more than 95% of genuine detections.This fact is illustrated and discussed in more detail inSection 7.1.In summary, when running
Duchamp on a realis-tic data cube with real galaxies at a redshift of about0 .
03 and genuine noise extracted from observationaldata taken with the WSRT, the software performs asexpected with completeness levels ranging in betweenthose achieved for the compact and extended modelsources discussed in the previous sections. This re-sult illustrates that the performance of
Duchamp , aswith any source finder based on flux thresholding, willstrongly depend on the morphology and extent of thesources to be detected. Even with multi-scale waveletreconstruction,
Duchamp is more likely to uncovercompact sources than sources that are significantly ex-tended, either spatially or spectrally.At the same time, the performance of
Duchamp does not seem to be hampered by the fact that we aredealing with real telescope data and noise, as the com-pleteness and reliability levels reported in Table 6 aregenerally very similar to what we achieved with themodel sources discussed in the previous sections. Thisis presumably due to the excellent quality of the West-erbork data which do not contain any obvious artefactssuch as interference or residual continuum emission.
In general,
Duchamp does what it promises to do. Itis able to reliably detect sources down to low signal-to-noise ratios and accurately determine their posi-tion and radial velocity. These are the most funda-mental requirements for any source finder. Our testsalso demonstrated that by using and fine-tuning theoptions of ‘`a trous’ wavelet reconstruction and grow-ing of sources to lower flux levels the performance of
Duchamp can be greatly enhanced.
The reliability figures reported throughout this pa-per have all been “raw” reliabilities, i.e. reliabilities asachieved by
Duchamp prior to any filtering of the out-put source catalogue. The user would normally wishto substantially improve these through appropriate fil-tering of the source catalogue based on the source pa-rameters as measured by
Duchamp .The left-hand panel of Figure 21 shows the mea-sured integrated flux plotted against measured linewidth for all genuine (black data points) and false(red data points) detections found by
Duchamp inthe point source models discussed in Section 4. Itis obvious that genuine and false detections occupylargely disjunct regions of F int – w parameter space,with false detections generally occurring near the lowend of the integrated flux spectrum. Similar plots canbe generated for other combinations of source param-eters, but F int and w usually provide the best dis-tinction between genuine and false detections.The easiest way to improve the reliability of Du-champ ’s source finding results is to simply apply a cutin F int to exclude most false detections while retainingmost of the genuine sources. In our example, apply-ing a cut at F int = 40 mJy km s − will discard 97 . .
9% of all genuine sources, thereby increasingthe overall reliability from 77 .
1% to 99 .
2% while onlymoderately decreasing the overall completeness from83 .
0% to 80 . ww.publish.csiro.au/journals/pasa C o m p l e t e n e ss / r e li a b ilit y ( % )
1D − 41D − 33D − 4 σσσ
Figure 19: Completeness (filled circles, solid lines)and reliability (open circles, dashed lines) of
Duchamp on the WSRT model cube with WHISPgalaxies for different flux thresholds and input pa-rameters. The colours, as shown in the legend,distinguish the different wavelet reconstructionmodes (one-dimensional versus three-dimensional)and wavelet reconstruction thresholds (3 σ versus4 σ ) used in the tests (see Table 6 for details). Thenumbers alongside the data points refer to the cor-responding runs as listed in Table 6. shifted WHISP galaxies, as plotted in the right-handpanel of Figure 21. Again, applying a simple fluxthreshold of 0 . − will improve reliability from10% to 84%, while only moderately decreasing com-pleteness from 59% to 50%. The method is not quiteas successful as for the point sources, as we are nowdealing with real galaxies and real noise with interfer-ence and artefacts, but nevertheless a significant im-provement in reliability can be achieved without anysevere impact on the number of genuine detections.This simple example illustrates that the “raw” reli-ability figures quoted throughout this paper should notbe considered as the final numbers. Reliability can begreatly improved through very basic filtering in param-eter space of the Duchamp output catalogue. In prin-ciple, this applies to the output of almost any sourcefinder. Alternatively, instead of removing sources fromthe output catalogue, it may be desirable to calculatea reliability number for each catalogue entry based onthe source’s location in parameter space and leave itto the catalogue’s users to decide as part of their sci-entific analysis at which reliability level they wish tomake the cut.
When it comes to source parametrisation, the mea-surements provided by
Duchamp are affected by sev-eral systematic errors. These systematic errors are notdue to errors in the software itself, but a consequenceof the the presence of noise in the data as well as themethods and algorithms used for measuring source pa-rameters. C o m p l e t e n e ss ( % ) Run 11Run 10Run 9Run 7Run 6Run 2
Figure 20: Completeness as a function of in-tegrated flux for selected runs (see legend) of
Duchamp on the WSRT model cube with WHISPgalaxies. The choice of colours is the same as inFigure 19.
Spectral line widths determined by
Duchamp aregenerally very accurate and not much affected by noise-induced, systematic errors as far as the w parameteris concerned. The two other line width parameterscalculated by Duchamp , w and w vel , appear to besystematically too large over a wide range of signal-to-noise ratios and should not be used unless explicitlyrequired in special, well-defined circumstances.Peak fluxes, as reported by Duchamp , are in gen-eral slightly too large for bright sources and signifi-cantly too large (on a relative scale) for faint sources.This is due to the fact that
Duchamp determines thepeak flux by simply selecting the value of the brightestpixel encountered. This method introduces a bias to-wards positive noise peaks sitting on top of the bright-est region of a source, and hence, in the presence ofnoise, peak fluxes measured by
Duchamp will be sys-tematically too high.Integrated fluxes determined by
Duchamp are sig-nificantly and systematically too small, in particularfor faint sources. This is likely caused by the fact that
Duchamp simply sums over the flux of discrete ele-ments above a given threshold to determine the in-tegrated flux, thereby missing some of the flux fromelements below the flux threshold. Hence, the rawintegrated flux measurements currently provided by
Duchamp are not useful and need to be corrected tocompensate for the systematic offset. This issue is par-ticularly sensitive as many scientific projects, includingthe ASKAP survey science projects WALLABY andDINGO (Meyer 2009), rely on accurate flux measure-ments, for example for determining the H i mass func-tion of galaxies.Finally, a particular problem in the case of galax-ies is that under certain circumstances galaxies eitherget broken up into multiple detections or only one half Deep Investigation of Neutral Gas Origins ; princi-pal investigator: Martin Meyer; public website: http://internal.physics.uwa.edu.au/ ∼ mmeyer/dingo/ Publications of the Astronomical Society of Australia M ea s u r e d i n t e g r a t e d f l ux ( J y k m / s ) M ea s u r e d i n t e g r a t e d f l ux ( J y k m / s ) WHISP galaxiesPoint sources
Figure 21: Measured integrated flux, F int , versus measured line width, w , of all genuine (black) andfalse (red) detections made by Duchamp in the point source models with Gaussian line profiles (left) andthe test cube with artificially redshifted WHISP galaxies (right). The dashed, black lines indicate the fluxlevels of 0 .
04 and 0 . − used to filter false detections. of a galaxy is detected. This problem mainly affectsfaint, edge-on galaxies with broad spectral profiles thatare partly hidden in the noise and results in systematicerrors in the measurements of essentially all source pa-rameters, including basic parameters such as positionand radial velocity. Such cases of multiple or partialdetections must be identified and treated separately toprevent biases in any scientific analysis based on thesource finding results. In this paper we present and discuss the results of ba-sic, three-dimensional source finding tests with
Du-champ , the standard source finder for the AustralianSKA Pathfinder, using different sets of unresolved andextended H i model sources as well as a data set of realgalaxies and noise obtained from H i observations withthe WSRT.Overall, Duchamp appears to be a successful, gen-eral-purpose source finder capable of reliably detect-ing sources down to low signal-to-noise ratios and ac-curately determining their position and velocity. Inthe case of point sources with simple Gaussian spec-tral lines we achieve a completeness of about 50% at apeak signal-to-noise ratio of 3 and an integrated fluxlevel of about 0 . − . The latter correspondsto an H i mass sensitivity of about 2 × M ⊙ at adistance of 100 Mpc which is slightly better than whatthe WALLABY project is expected to achieve for realgalaxies (Koribalski & Staveley-Smith 2009). The sit-uation is less ideal for extended sources with double-horn profiles. In this case we achieve 50% completenessat an integrated flux level of about 2 . − forthe model galaxies and 0 . − for the WHISPgalaxies. The latter is equivalent to an H i mass sensi-tivity of about 1 . × M ⊙ at a distance of 100 Mpc,illustrating that the performance of Duchamp , as wellas any other source finder, will strongly depend on source morphology. However, these figures may wellbe improved by carefully optimising the various inputparameters offered by
Duchamp .In its current state
Duchamp is not particularlysuccessful in parametrising sources in the presence ofnoise in the data cube, and other, external algorithmsfor source parametrisation should be considered in-stead. It appears, however, that most, if not all, para-metrisation issues are due to intrinsic limitations in theimplemented algorithms themselves and not due to er-rors in their implementation, suggesting that most ofthe problems can in principle be solved by implement-ing more sophisticated parametrisation algorithms in
Duchamp . Alternatively, corrections would have to beapplied to all parameters derived by
Duchamp to com-pensate for systematic errors. Such corrections, how-ever, would have to be highly specialised and tailoredto the particular survey and source type concerned.
Acknowledgments
We wish to thank Matthew Whiting for creating the
Duchamp source finder and for having had an openear for our comments and suggestions regarding im-provement of the software.
References
Barnes, D. G., et al. 2001, MNRAS, 322, 486DeBoer, D. R., et al. 2009, IEEE Proceedings, 97, 1507Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio,T. J. L. 2009, IEEE Proceedings, 97, 1482Giovanelli, R. 2005, AJ, 130, 2598Jonas, J. L. 2009, IEEE Proceedings, 97, 1522 ww.publish.csiro.au/journals/pasa this issue Saintonge, A. 2007, AJ, 133, 2087Sault, R. J., Teuben, P. J., & Wright M. C. H.1995, in: Astronomical Data Analysis Softwareand Systems IV, eds. R. Shaw, H. E. Payne, &J. J. E. Hayes, ASP Conf. Ser., 77, 433Serra, P., Jurek, R., & Fl¨oer, L. 2011, this issue
Starck, J.-L. & Murtagh, F. Astronomical Image andData Analysis, Springer, 2002Swaters, R. A., van Albada, T. S., van der Hulst, J. M.,& Sancisi, R. 2002, A&A, 390, 829Westmeier, T. & Johnston, S. 2010, in: InternationalSKA Forum 2010 Science Meeting, Proceedings ofScience, PoS(ISKAF2010)056Whiting, M. T. 2011a, this issue
Whiting, M. T. 2011b,