[303] | 1 | % ----------------------------------------------------------------------- |
---|
| 2 | % executionFlow.tex: Section detailing each of the main algorithms |
---|
| 3 | % used by Duchamp. |
---|
| 4 | % ----------------------------------------------------------------------- |
---|
| 5 | % Copyright (C) 2006, Matthew Whiting, ATNF |
---|
| 6 | % |
---|
| 7 | % This program is free software; you can redistribute it and/or modify it |
---|
| 8 | % under the terms of the GNU General Public License as published by the |
---|
| 9 | % Free Software Foundation; either version 2 of the License, or (at your |
---|
| 10 | % option) any later version. |
---|
| 11 | % |
---|
| 12 | % Duchamp is distributed in the hope that it will be useful, but WITHOUT |
---|
| 13 | % ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
---|
| 14 | % FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
| 15 | % for more details. |
---|
| 16 | % |
---|
| 17 | % You should have received a copy of the GNU General Public License |
---|
| 18 | % along with Duchamp; if not, write to the Free Software Foundation, |
---|
| 19 | % Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA |
---|
| 20 | % |
---|
| 21 | % Correspondence concerning Duchamp may be directed to: |
---|
| 22 | % Internet email: Matthew.Whiting [at] atnf.csiro.au |
---|
| 23 | % Postal address: Dr. Matthew Whiting |
---|
| 24 | % Australia Telescope National Facility, CSIRO |
---|
| 25 | % PO Box 76 |
---|
| 26 | % Epping NSW 1710 |
---|
| 27 | % AUSTRALIA |
---|
| 28 | % ----------------------------------------------------------------------- |
---|
[258] | 29 | \secA{What \duchamp is doing} |
---|
[158] | 30 | \label{sec-flow} |
---|
| 31 | |
---|
[265] | 32 | Each of the steps that \duchamp goes through in the course of its |
---|
| 33 | execution are discussed here in more detail. This should provide |
---|
| 34 | enough background information to fully understand what \duchamp is |
---|
| 35 | doing and what all the output information is. For those interested in |
---|
| 36 | the programming side of things, \duchamp is written in C/C++ and makes |
---|
| 37 | use of the \textsc{cfitsio}, \textsc{wcslib} and \textsc{pgplot} |
---|
[158] | 38 | libraries. |
---|
| 39 | |
---|
| 40 | \secB{Image input} |
---|
| 41 | \label{sec-input} |
---|
| 42 | |
---|
[162] | 43 | The cube is read in using basic \textsc{cfitsio} commands, and stored |
---|
| 44 | as an array in a special C++ class. This class keeps track of the list |
---|
| 45 | of detected objects, as well as any reconstructed arrays that are made |
---|
| 46 | (see \S\ref{sec-recon}). The World Coordinate System |
---|
| 47 | (WCS)\footnote{This is the information necessary for translating the |
---|
| 48 | pixel locations to quantities such as position on the sky, frequency, |
---|
| 49 | velocity, and so on.} information for the cube is also obtained from |
---|
| 50 | the FITS header by \textsc{wcslib} functions \citep{greisen02, |
---|
| 51 | calabretta02}, and this information, in the form of a \texttt{wcsprm} |
---|
| 52 | structure, is also stored in the same class. |
---|
[158] | 53 | |
---|
[231] | 54 | A sub-section of a cube can be requested by defining the subsection |
---|
| 55 | with the \texttt{subsection} parameter and setting |
---|
[298] | 56 | \texttt{flagSubsection = true} -- this can be a good idea if the cube |
---|
[231] | 57 | has very noisy edges, which may produce many spurious detections. |
---|
| 58 | |
---|
| 59 | There are two ways of specifying the \texttt{subsection} string. The |
---|
| 60 | first is the generalised form |
---|
| 61 | \texttt{[x1:x2:dx,y1:y2:dy,z1:z2:dz,...]}, as used by the |
---|
| 62 | \textsc{cfitsio} library. This has one set of colon-separated numbers |
---|
| 63 | for each axis in the FITS file. In this manner, the x-coordinates run |
---|
[158] | 64 | from \texttt{x1} to \texttt{x2} (inclusive), with steps of |
---|
[231] | 65 | \texttt{dx}. The step value can be omitted, so a subsection of the |
---|
[258] | 66 | form \texttt{[2:50,2:50,10:1000]} is still valid. In fact, \duchamp |
---|
[231] | 67 | does not make use of any step value present in the subsection string, |
---|
| 68 | and any that are present are removed before the file is opened. |
---|
[158] | 69 | |
---|
[231] | 70 | If the entire range of a coordinate is required, one can replace the |
---|
| 71 | range with a single asterisk, \eg \texttt{[2:50,2:50,*]}. Thus, the |
---|
[298] | 72 | subsection string \texttt{[*,*,*]} is simply the entire cube. Note |
---|
| 73 | that the pixel ranges for each axis start at 1, so the full pixel |
---|
| 74 | range of a 100-pixel axis would be expressed as 1:100. A complete |
---|
| 75 | description of this section syntax can be found at the |
---|
[231] | 76 | \textsc{fitsio} web site% |
---|
[158] | 77 | \footnote{% |
---|
| 78 | \href% |
---|
[223] | 79 | {http://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c\_user/node91.html}% |
---|
| 80 | {http://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c\_user/node91.html}}. |
---|
[158] | 81 | |
---|
[231] | 82 | |
---|
| 83 | Making full use of the subsection requires knowledge of the size of |
---|
| 84 | each of the dimensions. If one wants to, for instance, trim a certain |
---|
| 85 | number of pixels off the edges of the cube, without examining the cube |
---|
| 86 | to obtain the actual size, one can use the second form of the |
---|
| 87 | subsection string. This just gives a number for each axis, \eg |
---|
| 88 | \texttt{[5,5,5]} (which would trim 5 pixels from the start \emph{and} |
---|
| 89 | end of each axis). |
---|
| 90 | |
---|
[298] | 91 | All types of subsections can be combined \eg \texttt{[5,2:98,*]}. |
---|
[231] | 92 | |
---|
[447] | 93 | Typically, the units of pixel brightness are given by the FITS file's |
---|
| 94 | BUNIT keyword. However, this may often be unwieldy (for instance, the |
---|
| 95 | units are Jy/beam, but the values are around a few mJy/beam). It is |
---|
| 96 | therefore possible to nominate new units, to which the pixel values |
---|
| 97 | will be converted, by using the \texttt{newFluxUnits} input |
---|
| 98 | parameter. The units must be directly translatable from the existing |
---|
| 99 | ones -- for instance, if BUNIT is Jy/beam, you cannot specify mJy, it |
---|
| 100 | must be mJy/beam. If an incompatible unit is given, the BUNIT value is |
---|
| 101 | used instead. |
---|
[231] | 102 | |
---|
[158] | 103 | \secB{Image modification} |
---|
| 104 | \label{sec-modify} |
---|
| 105 | |
---|
| 106 | Several modifications to the cube can be made that improve the |
---|
[258] | 107 | execution and efficiency of \duchamp (their use is optional, governed |
---|
[158] | 108 | by the relevant flags in the parameter file). |
---|
| 109 | |
---|
| 110 | \secC{BLANK pixel removal} |
---|
[285] | 111 | \label{sec-blank} |
---|
[158] | 112 | |
---|
[162] | 113 | If the imaged area of a cube is non-rectangular (see the example in |
---|
[285] | 114 | Fig.~\ref{fig-moment}, a cube from the HIPASS survey), BLANK pixels |
---|
| 115 | are used to pad it out to a rectangular shape. The value of these |
---|
| 116 | pixels is given by the FITS header keywords BLANK, BSCALE and |
---|
| 117 | BZERO. While these pixels make the image a nice shape, they will take |
---|
| 118 | up unnecessary space in memory, and so to potentially speed up the |
---|
| 119 | processing we can trim them from the edge. This is done when the |
---|
[298] | 120 | parameter \texttt{flagTrim = true}. If the above keywords are not |
---|
[285] | 121 | present, the trimming will not be done (in this case, a similar effect |
---|
| 122 | can be accomplished, if one knows where the ``blank'' pixels are, by |
---|
| 123 | using the subsection option). |
---|
[158] | 124 | |
---|
[285] | 125 | The amount of trimming is recorded, and these pixels are added back in |
---|
| 126 | once the source-detection is completed (so that quoted pixel positions |
---|
| 127 | are applicable to the original cube). Rows and columns are trimmed one |
---|
| 128 | at a time until the first non-BLANK pixel is reached, so that the |
---|
| 129 | image remains rectangular. In practice, this means that there will be |
---|
| 130 | some BLANK pixels left in the trimmed image (if the non-BLANK region |
---|
| 131 | is non-rectangular). However, these are ignored in all further |
---|
| 132 | calculations done on the cube. |
---|
[158] | 133 | |
---|
| 134 | \secC{Baseline removal} |
---|
| 135 | |
---|
| 136 | Second, the user may request the removal of baselines from the |
---|
| 137 | spectra, via the parameter \texttt{flagBaseline}. This may be |
---|
| 138 | necessary if there is a strong baseline ripple present, which can |
---|
| 139 | result in spurious detections at the high points of the ripple. The |
---|
| 140 | baseline is calculated from a wavelet reconstruction procedure (see |
---|
| 141 | \S\ref{sec-recon}) that keeps only the two largest scales. This is |
---|
| 142 | done separately for each spatial pixel (\ie for each spectrum in the |
---|
| 143 | cube), and the baselines are stored and added back in before any |
---|
| 144 | output is done. In this way the quoted fluxes and displayed spectra |
---|
| 145 | are as one would see from the input cube itself -- even though the |
---|
| 146 | detection (and reconstruction if applicable) is done on the |
---|
| 147 | baseline-removed cube. |
---|
| 148 | |
---|
| 149 | The presence of very strong signals (for instance, masers at several |
---|
[162] | 150 | hundred Jy) could affect the determination of the baseline, and would |
---|
| 151 | lead to a large dip centred on the signal in the baseline-subtracted |
---|
[158] | 152 | spectrum. To prevent this, the signal is trimmed prior to the |
---|
| 153 | reconstruction process at some standard threshold (at $8\sigma$ above |
---|
| 154 | the mean). The baseline determined should thus be representative of |
---|
| 155 | the true, signal-free baseline. Note that this trimming is only a |
---|
| 156 | temporary measure which does not affect the source-detection. |
---|
| 157 | |
---|
| 158 | \secC{Ignoring bright Milky Way emission} |
---|
[386] | 159 | \label{sec-MW} |
---|
[158] | 160 | |
---|
| 161 | Finally, a single set of contiguous channels can be ignored -- these |
---|
| 162 | may exhibit very strong emission, such as that from the Milky Way as |
---|
[258] | 163 | seen in extragalactic \hi cubes (hence the references to ``Milky |
---|
[158] | 164 | Way'' in relation to this task -- apologies to Galactic |
---|
| 165 | astronomers!). Such dominant channels will produce many detections |
---|
| 166 | that are unnecessary, uninteresting (if one is interested in |
---|
| 167 | extragalactic \hi) and large (in size and hence in memory usage), and |
---|
| 168 | so will slow the program down and detract from the interesting |
---|
| 169 | detections. |
---|
| 170 | |
---|
| 171 | The use of this feature is controlled by the \texttt{flagMW} |
---|
| 172 | parameter, and the exact channels concerned are able to be set by the |
---|
| 173 | user (using \texttt{maxMW} and \texttt{minMW} -- these give an |
---|
| 174 | inclusive range of channels). When employed, these channels are |
---|
| 175 | ignored for the searching, and the scaling of the spectral output (see |
---|
| 176 | Fig.~\ref{fig-spect}) will not take them into account. They will be |
---|
| 177 | present in the reconstructed array, however, and so will be included |
---|
| 178 | in the saved FITS file (see \S\ref{sec-reconIO}). When the final |
---|
| 179 | spectra are plotted, the range of channels covered by these parameters |
---|
[673] | 180 | is indicated by a green hashed box. Note that these channels refer to |
---|
| 181 | channel numbers in the full cube, before any subsection is applied. |
---|
[158] | 182 | |
---|
| 183 | \secB{Image reconstruction} |
---|
| 184 | \label{sec-recon} |
---|
| 185 | |
---|
[258] | 186 | The user can direct \duchamp to reconstruct the data cube using the |
---|
| 187 | \atrous wavelet procedure. A good description of the procedure can be |
---|
[158] | 188 | found in \citet{starck02:book}. The reconstruction is an effective way |
---|
| 189 | of removing a lot of the noise in the image, allowing one to search |
---|
| 190 | reliably to fainter levels, and reducing the number of spurious |
---|
| 191 | detections. This is an optional step, but one that greatly enhances |
---|
| 192 | the source-detection process, with the payoff that it can be |
---|
| 193 | relatively time- and memory-intensive. |
---|
| 194 | |
---|
| 195 | \secC{Algorithm} |
---|
| 196 | |
---|
[258] | 197 | The steps in the \atrous reconstruction are as follows: |
---|
[158] | 198 | \begin{enumerate} |
---|
[162] | 199 | \item The reconstructed array is set to 0 everywhere. |
---|
[158] | 200 | \item The input array is discretely convolved with a given filter |
---|
| 201 | function. This is determined from the parameter file via the |
---|
| 202 | \texttt{filterCode} parameter -- see Appendix~\ref{app-param} for |
---|
| 203 | details on the filters available. |
---|
| 204 | \item The wavelet coefficients are calculated by taking the difference |
---|
| 205 | between the convolved array and the input array. |
---|
| 206 | \item If the wavelet coefficients at a given point are above the |
---|
| 207 | requested threshold (given by \texttt{snrRecon} as the number of |
---|
| 208 | $\sigma$ above the mean and adjusted to the current scale -- see |
---|
| 209 | Appendix~\ref{app-scaling}), add these to the reconstructed array. |
---|
[285] | 210 | \item The separation between the filter coefficients is doubled. (Note |
---|
| 211 | that this step provides the name of the procedure\footnote{\atrous |
---|
| 212 | means ``with holes'' in French.}, as gaps or holes are created in |
---|
| 213 | the filter coverage.) |
---|
[158] | 214 | \item The procedure is repeated from step 2, using the convolved array |
---|
| 215 | as the input array. |
---|
| 216 | \item Continue until the required maximum number of scales is reached. |
---|
| 217 | \item Add the final smoothed (\ie convolved) array to the |
---|
| 218 | reconstructed array. This provides the ``DC offset'', as each of the |
---|
| 219 | wavelet coefficient arrays will have zero mean. |
---|
| 220 | \end{enumerate} |
---|
| 221 | |
---|
[364] | 222 | The range of scales at which the selection of wavelet coefficients is |
---|
| 223 | made is governed by the \texttt{scaleMin} and \texttt{scaleMax} |
---|
| 224 | parameters. The minimum scale used is given by \texttt{scaleMin}, |
---|
| 225 | where the default value is 1 (the first scale). This parameter is |
---|
| 226 | useful if you want to ignore the highest-frequency features |
---|
| 227 | (e.g. high-frequency noise that might be present). Normally the |
---|
| 228 | maximum scale is calculated from the size of the input array, but it |
---|
| 229 | can be specified by using \texttt{scaleMax}. A value $\le0$ will |
---|
| 230 | result in the use of the calculated value, as will a value of |
---|
| 231 | \texttt{scaleMax} greater than the calculated value. Use of these two |
---|
| 232 | parameters can allow searching for features of a particular scale |
---|
| 233 | size, for instance searching for narrow absorption features. |
---|
| 234 | |
---|
[158] | 235 | The reconstruction has at least two iterations. The first iteration |
---|
| 236 | makes a first pass at the wavelet reconstruction (the process outlined |
---|
[162] | 237 | in the 8 stages above), but the residual array will likely have some |
---|
| 238 | structure still in it, so the wavelet filtering is done on the |
---|
[158] | 239 | residual, and any significant wavelet terms are added to the final |
---|
[162] | 240 | reconstruction. This step is repeated until the change in the measured |
---|
| 241 | standard deviation of the background (see note below on the evaluation |
---|
| 242 | of this quantity) is less than some fiducial amount. |
---|
[158] | 243 | |
---|
[258] | 244 | It is important to note that the \atrous decomposition is an example |
---|
[158] | 245 | of a ``redundant'' transformation. If no thresholding is performed, |
---|
| 246 | the sum of all the wavelet coefficient arrays and the final smoothed |
---|
| 247 | array is identical to the input array. The thresholding thus removes |
---|
| 248 | only the unwanted structure in the array. |
---|
| 249 | |
---|
| 250 | Note that any BLANK pixels that are still in the cube will not be |
---|
| 251 | altered by the reconstruction -- they will be left as BLANK so that |
---|
| 252 | the shape of the valid part of the cube is preserved. |
---|
| 253 | |
---|
| 254 | \secC{Note on Statistics} |
---|
| 255 | |
---|
| 256 | The correct calculation of the reconstructed array needs good |
---|
[386] | 257 | estimators of the underlying mean and standard deviation (or rms) of |
---|
| 258 | the background noise distribution. The methods used to estimate these |
---|
| 259 | quantities are detailed in \S\ref{sec-stats} -- the default behaviour |
---|
| 260 | is to use robust estimators, to avoid biasing due to bright pixels. |
---|
[158] | 261 | |
---|
[386] | 262 | %These statistics are estimated using |
---|
| 263 | %robust methods, to avoid corruption by strong outlying points. The |
---|
| 264 | %mean of the distribution is actually estimated by the median, while |
---|
| 265 | %the median absolute deviation from the median (MADFM) is calculated |
---|
| 266 | %and corrected assuming Gaussianity to estimate the underlying standard |
---|
| 267 | %deviation $\sigma$. The Gaussianity (or Normality) assumption is |
---|
| 268 | %critical, as the MADFM does not give the same value as the usual rms |
---|
| 269 | %or standard deviation value -- for a Normal distribution |
---|
| 270 | %$N(\mu,\sigma)$ we find MADFM$=0.6744888\sigma$, but this will change |
---|
| 271 | %for different distributions. Since this ratio is corrected for, the |
---|
| 272 | %user need only think in the usual multiples of the rms when setting |
---|
| 273 | %\texttt{snrRecon}. See Appendix~\ref{app-madfm} for a derivation of |
---|
| 274 | %this value. |
---|
| 275 | |
---|
[265] | 276 | When thresholding the different wavelet scales, the value of the rms |
---|
[158] | 277 | as measured from the wavelet array needs to be scaled to account for |
---|
| 278 | the increased amount of correlation between neighbouring pixels (due |
---|
| 279 | to the convolution). See Appendix~\ref{app-scaling} for details on |
---|
| 280 | this scaling. |
---|
| 281 | |
---|
| 282 | \secC{User control of reconstruction parameters} |
---|
| 283 | |
---|
| 284 | The most important parameter for the user to select in relation to the |
---|
| 285 | reconstruction is the threshold for each wavelet array. This is set |
---|
| 286 | using the \texttt{snrRecon} parameter, and is given as a multiple of |
---|
| 287 | the rms (estimated by the MADFM) above the mean (which for the wavelet |
---|
| 288 | arrays should be approximately zero). There are several other |
---|
| 289 | parameters that can be altered as well that affect the outcome of the |
---|
| 290 | reconstruction. |
---|
| 291 | |
---|
| 292 | By default, the cube is reconstructed in three dimensions, using a |
---|
| 293 | 3-dimensional filter and 3-dimensional convolution. This can be |
---|
| 294 | altered, however, using the parameter \texttt{reconDim}. If set to 1, |
---|
| 295 | this means the cube is reconstructed by considering each spectrum |
---|
| 296 | separately, whereas \texttt{reconDim=2} will mean the cube is |
---|
| 297 | reconstructed by doing each channel map separately. The merits of |
---|
| 298 | these choices are discussed in \S\ref{sec-notes}, but it should be |
---|
| 299 | noted that a 2-dimensional reconstruction can be susceptible to edge |
---|
[162] | 300 | effects if the spatial shape of the pixel array is not rectangular. |
---|
[158] | 301 | |
---|
| 302 | The user can also select the minimum scale to be used in the |
---|
| 303 | reconstruction. The first scale exhibits the highest frequency |
---|
| 304 | variations, and so ignoring this one can sometimes be beneficial in |
---|
| 305 | removing excess noise. The default is to use all scales |
---|
| 306 | (\texttt{minscale = 1}). |
---|
| 307 | |
---|
| 308 | Finally, the filter that is used for the convolution can be selected |
---|
| 309 | by using \texttt{filterCode} and the relevant code number -- the |
---|
| 310 | choices are listed in Appendix~\ref{app-param}. A larger filter will |
---|
| 311 | give a better reconstruction, but take longer and use more memory when |
---|
| 312 | executing. When multi-dimensional reconstruction is selected, this |
---|
| 313 | filter is used to construct a 2- or 3-dimensional equivalent. |
---|
| 314 | |
---|
[208] | 315 | \secB{Smoothing the cube} |
---|
| 316 | \label{sec-smoothing} |
---|
| 317 | |
---|
[275] | 318 | An alternative to doing the wavelet reconstruction is to smooth the |
---|
| 319 | cube. This technique can be useful in reducing the noise level |
---|
| 320 | slightly (at the cost of making neighbouring pixels correlated and |
---|
| 321 | blurring any signal present), and is particularly well suited to the |
---|
| 322 | case where a particular signal size (\ie a certain channel width or |
---|
| 323 | spatial size) is believed to be present in the data. |
---|
[208] | 324 | |
---|
[275] | 325 | There are two alternative methods that can be used: spectral |
---|
| 326 | smoothing, using the Hanning filter; or spatial smoothing, using a 2D |
---|
| 327 | Gaussian kernel. These alternatives are outlined below. To utilise the |
---|
| 328 | smoothing option, set the parameter \texttt{flagSmooth=true} and set |
---|
| 329 | \texttt{smoothType} to either \texttt{spectral} or \texttt{spatial}. |
---|
[208] | 330 | |
---|
[275] | 331 | \secC{Spectral smoothing} |
---|
| 332 | |
---|
[298] | 333 | When \texttt{smoothType = spectral} is selected, the cube is smoothed |
---|
[275] | 334 | only in the spectral domain. Each spectrum is independently smoothed |
---|
| 335 | by a Hanning filter, and then put back together to form the smoothed |
---|
| 336 | cube, which is then used by the searching algorithm (see below). Note |
---|
| 337 | that in the case of both the reconstruction and the smoothing options |
---|
| 338 | being requested, the reconstruction will take precedence and the |
---|
| 339 | smoothing will \emph{not} be done. |
---|
| 340 | |
---|
[208] | 341 | There is only one parameter necessary to define the degree of |
---|
| 342 | smoothing -- the Hanning width $a$ (given by the user parameter |
---|
[231] | 343 | \texttt{hanningWidth}). The coefficients $c(x)$ of the Hanning filter |
---|
| 344 | are defined by |
---|
[208] | 345 | \[ |
---|
[231] | 346 | c(x) = |
---|
| 347 | \begin{cases} |
---|
[645] | 348 | \frac{1}{2}\left(1+\cos(\frac{\pi x}{a})\right) &|x| < (a+1)/2\\ |
---|
| 349 | 0 &|x| \geq (a+1)/2. |
---|
[285] | 350 | \end{cases},\ a,x \in \mathbb{Z} |
---|
[208] | 351 | \] |
---|
[277] | 352 | Note that the width specified must be an |
---|
[232] | 353 | odd integer (if the parameter provided is even, it is incremented by |
---|
| 354 | one). |
---|
[208] | 355 | |
---|
[275] | 356 | \secC{Spatial smoothing} |
---|
[208] | 357 | |
---|
[298] | 358 | When \texttt{smoothType = spatial} is selected, the cube is smoothed |
---|
[275] | 359 | only in the spatial domain. Each channel map is independently smoothed |
---|
[285] | 360 | by a two-dimensional Gaussian kernel, put back together to form the |
---|
| 361 | smoothed cube, and used in the searching algorithm (see below). Again, |
---|
| 362 | reconstruction is always done by preference if both techniques are |
---|
| 363 | requested. |
---|
[275] | 364 | |
---|
| 365 | The two-dimensional Gaussian has three parameters to define it, |
---|
[285] | 366 | governed by the elliptical cross-sectional shape of the Gaussian |
---|
[275] | 367 | function: the FWHM (full-width at half-maximum) of the major and minor |
---|
| 368 | axes, and the position angle of the major axis. These are given by the |
---|
[298] | 369 | user parameters \texttt{kernMaj, kernMin} \& \texttt{kernPA}. If a |
---|
| 370 | circular Gaussian is required, the user need only provide the |
---|
| 371 | \texttt{kernMaj} parameter. The \texttt{kernMin} parameter will then |
---|
| 372 | be set to the same value, and \texttt{kernPA} to zero. If we define |
---|
| 373 | these parameters as $a,b,\theta$ respectively, we can define the |
---|
| 374 | kernel by the function |
---|
[275] | 375 | \[ |
---|
[277] | 376 | k(x,y) = \exp\left[-0.5 \left(\frac{X^2}{\sigma_X^2} + |
---|
| 377 | \frac{Y^2}{\sigma_Y^2} \right) \right] |
---|
[275] | 378 | \] |
---|
| 379 | where $(x,y)$ are the offsets from the central pixel of the gaussian |
---|
| 380 | function, and |
---|
[277] | 381 | \begin{align*} |
---|
| 382 | X& = x\sin\theta - y\cos\theta& |
---|
| 383 | Y&= x\cos\theta + y\sin\theta\\ |
---|
| 384 | \sigma_X^2& = \frac{(a/2)^2}{2\ln2}& |
---|
| 385 | \sigma_Y^2& = \frac{(b/2)^2}{2\ln2}\\ |
---|
| 386 | \end{align*} |
---|
[275] | 387 | |
---|
[285] | 388 | \secB{Input/Output of reconstructed/smoothed arrays} |
---|
[277] | 389 | \label{sec-reconIO} |
---|
| 390 | |
---|
| 391 | The smoothing and reconstruction stages can be relatively |
---|
| 392 | time-consuming, particularly for large cubes and reconstructions in |
---|
| 393 | 3-D (or even spatial smoothing). To get around this, \duchamp provides |
---|
| 394 | a shortcut to allow users to perform multiple searches (\eg with |
---|
| 395 | different thresholds) on the same reconstruction/smoothing setup |
---|
| 396 | without re-doing the calculations each time. |
---|
| 397 | |
---|
| 398 | To save the reconstructed array as a FITS file, set |
---|
| 399 | \texttt{flagOutputRecon = true}. The file will be saved in the same |
---|
| 400 | directory as the input image, so the user needs to have write |
---|
| 401 | permissions for that directory. |
---|
| 402 | |
---|
[525] | 403 | The name of the file can given by the \texttt{fileOutputRecon} |
---|
| 404 | parameter, but this can be ignored and \duchamp will present a name |
---|
| 405 | based on the reconstruction parameters. The filename will be derived |
---|
| 406 | from the input filename, with extra information detailing the |
---|
| 407 | reconstruction that has been done. For example, suppose |
---|
| 408 | \texttt{image.fits} has been reconstructed using a 3-dimensional |
---|
| 409 | reconstruction with filter \#2, thresholded at $4\sigma$ using all |
---|
| 410 | scales. The output filename will then be |
---|
[277] | 411 | \texttt{image.RECON-3-2-4-1.fits} (\ie it uses the four parameters |
---|
| 412 | relevant for the \atrous reconstruction as listed in |
---|
| 413 | Appendix~\ref{app-param}). The new FITS file will also have these |
---|
| 414 | parameters as header keywords. If a subsection of the input image has |
---|
| 415 | been used (see \S\ref{sec-input}), the format of the output filename |
---|
| 416 | will be \texttt{image.sub.RECON-3-2-4-1.fits}, and the subsection that |
---|
| 417 | has been used is also stored in the FITS header. |
---|
| 418 | |
---|
| 419 | Likewise, the residual image, defined as the difference between the |
---|
| 420 | input and reconstructed arrays, can also be saved in the same manner |
---|
| 421 | by setting \texttt{flagOutputResid = true}. Its filename will be the |
---|
| 422 | same as above, with \texttt{RESID} replacing \texttt{RECON}. |
---|
| 423 | |
---|
| 424 | If a reconstructed image has been saved, it can be read in and used |
---|
| 425 | instead of redoing the reconstruction. To do so, the user should set |
---|
| 426 | the parameter \texttt{flagReconExists = true}. The user can indicate |
---|
| 427 | the name of the reconstructed FITS file using the \texttt{reconFile} |
---|
| 428 | parameter, or, if this is not specified, \duchamp searches for the |
---|
| 429 | file with the name as defined above. If the file is not found, the |
---|
| 430 | reconstruction is performed as normal. Note that to do this, the user |
---|
| 431 | needs to set \texttt{flagAtrous = true} (obviously, if this is |
---|
| 432 | \texttt{false}, the reconstruction is not needed). |
---|
| 433 | |
---|
[525] | 434 | To save the smoothed array, set \texttt{flagOutputSmooth = true}. As |
---|
| 435 | for the reconstructed/residual arrays, the |
---|
| 436 | name of the file can given by the \texttt{fileOutputSmooth} parameter, |
---|
| 437 | but this can be ignored and \duchamp will present a name based on the |
---|
| 438 | method of smoothing used. It will be either |
---|
| 439 | \texttt{image.SMOOTH-1D-a.fits}, where a is replaced by the Hanning |
---|
| 440 | width used, or \texttt{image.SMOOTH-2D-a-b-c.fits}, where the Gaussian |
---|
| 441 | kernel parameters are a,b,c. Similarly to the reconstruction case, a |
---|
| 442 | saved file can be read in by setting \texttt{flagSmoothExists = true} |
---|
| 443 | and either specifying a file to be read with the \texttt{smoothFile} |
---|
| 444 | parameter or relying on \duchamp to find the file with the name as |
---|
| 445 | given above. |
---|
[277] | 446 | |
---|
| 447 | |
---|
[158] | 448 | \secB{Searching the image} |
---|
| 449 | \label{sec-detection} |
---|
| 450 | |
---|
[277] | 451 | \secC{Technique} |
---|
| 452 | |
---|
[298] | 453 | The basic idea behind detection in \duchamp is to locate sets of |
---|
| 454 | contiguous voxels that lie above some threshold. No size or shape |
---|
| 455 | requirement is imposed upon the detections -- that is, \duchamp does |
---|
| 456 | not fit \eg a Gaussian profile to each source. All it does is find |
---|
| 457 | connected groups of bright voxels. |
---|
[258] | 458 | |
---|
[298] | 459 | One threshold is calculated for the entire cube, enabling calculation |
---|
| 460 | of signal-to-noise ratios for each source (see |
---|
| 461 | Section~\ref{sec-output} for details). The user can manually specify a |
---|
| 462 | value (using the parameter \texttt{threshold}) for the threshold, |
---|
[462] | 463 | which will override the calculated value. Note that this option |
---|
| 464 | overrides any settings of \texttt{snrCut} or FDR options (see below). |
---|
[298] | 465 | |
---|
[686] | 466 | The cube can be searched in one of two ways, governed by the input |
---|
| 467 | parameter \texttt{searchType}. If \texttt{searchType=spatial}, the |
---|
| 468 | cube is searched one channel map at a time, using the 2-dimensional |
---|
| 469 | raster-scanning algorithm of \citet{lutz80} that connects groups of |
---|
| 470 | neighbouring pixels. Such an algorithm cannot be applied directly to a |
---|
| 471 | 3-dimensional case, as it requires that objects are completely nested |
---|
| 472 | in a row (when scanning along a row, if an object finishes and other |
---|
| 473 | starts, you won't get back to the first until the second is completely |
---|
| 474 | finished for the row). Three-dimensional data does not have this |
---|
| 475 | property, hence the need to treat the data on a 2-dimensional basis at |
---|
| 476 | most. |
---|
[158] | 477 | |
---|
[686] | 478 | Alternatively, if \texttt{searchType=spectral}, the searching is done |
---|
| 479 | in one dimension on each individual spatial pixel's spectrum. This is |
---|
| 480 | a simpler search, but there are potentially many more of them. |
---|
| 481 | |
---|
[265] | 482 | Although there are parameters that govern the minimum number of pixels |
---|
[720] | 483 | in a spatial, spectral and total senses that an object must have |
---|
| 484 | (\texttt{minPix}, \texttt{minChannels} and \texttt{minVoxels} |
---|
| 485 | respectively), these criteria are not applied at this point - see |
---|
| 486 | \S\ref{sec-reject} for details. |
---|
[158] | 487 | |
---|
[258] | 488 | Finally, the search only looks for positive features. If one is |
---|
| 489 | interested instead in negative features (such as absorption lines), |
---|
| 490 | set the parameter \texttt{flagNegative = true}. This will invert the |
---|
| 491 | cube (\ie multiply all pixels by $-1$) prior to the search, and then |
---|
| 492 | re-invert the cube (and the fluxes of any detections) after searching |
---|
| 493 | is complete. All outputs are done in the same manner as normal, so |
---|
| 494 | that fluxes of detections will be negative. |
---|
[158] | 495 | |
---|
[258] | 496 | \secC{Calculating statistics} |
---|
[386] | 497 | \label{sec-stats} |
---|
[258] | 498 | |
---|
[386] | 499 | A crucial part of the detection process (as well as the wavelet |
---|
| 500 | reconstruction: \S\ref{sec-recon}) is estimating the statistics that |
---|
| 501 | define the detection threshold. To determine a threshold, we need to |
---|
| 502 | estimate from the data two parameters: the middle of the noise |
---|
[277] | 503 | distribution (the ``noise level''), and the width of the distribution |
---|
[386] | 504 | (the ``noise spread''). The noise level is estimated by either the |
---|
| 505 | mean or the median, and the noise spread by the rms (or the standard |
---|
| 506 | deviation) or the median absolute deviation from the median |
---|
| 507 | (MADFM). The median and MADFM are robust statistics, in that they are |
---|
| 508 | not biased by the presence of a few pixels much brighter than the |
---|
| 509 | noise. |
---|
[258] | 510 | |
---|
[386] | 511 | All four statistics are calculated automatically, but the choice of |
---|
| 512 | parameters that will be used is governed by the input parameter |
---|
| 513 | \texttt{flagRobustStats}. This has the default value \texttt{true}, |
---|
| 514 | meaning the underlying mean of the noise distribution is estimated by |
---|
| 515 | the median, and the underlying standard deviation is estimated by the |
---|
| 516 | MADFM. In the latter case, the value is corrected, under the |
---|
| 517 | assumption that the underlying distribution is Normal (Gaussian), by |
---|
| 518 | dividing by 0.6744888 -- see Appendix~\ref{app-madfm} for details. If |
---|
| 519 | \texttt{flagRobustStats=false}, the mean and rms are used instead. |
---|
| 520 | |
---|
[277] | 521 | The choice of pixels to be used depend on the analysis method. If the |
---|
| 522 | wavelet reconstruction has been done, the residuals (defined |
---|
| 523 | in the sense of original $-$ reconstruction) are used to estimate the |
---|
| 524 | noise spread of the cube, since the reconstruction should pick out |
---|
| 525 | all significant structure. The noise level (the middle of the |
---|
| 526 | distribution) is taken from the original array. |
---|
| 527 | |
---|
| 528 | If smoothing of the cube has been done instead, all noise parameters |
---|
| 529 | are measured from the smoothed array, and detections are made with |
---|
| 530 | these parameters. When the signal-to-noise level is quoted for each |
---|
| 531 | detection (see \S\ref{sec-output}), the noise parameters of the |
---|
| 532 | original array are used, since the smoothing process correlates |
---|
| 533 | neighbouring pixels, reducing the noise level. |
---|
| 534 | |
---|
| 535 | If neither reconstruction nor smoothing has been done, then the |
---|
| 536 | statistics are calculated from the original, input array. |
---|
| 537 | |
---|
[258] | 538 | The parameters that are estimated should be representative of the |
---|
| 539 | noise in the cube. For the case of small objects embedded in many |
---|
| 540 | noise pixels (\eg the case of \hi surveys), using the full cube will |
---|
| 541 | provide good estimators. It is possible, however, to use only a |
---|
| 542 | subsection of the cube by setting the parameter \texttt{flagStatSec = |
---|
[386] | 543 | true} and providing the desired subsection to the \texttt{StatSec} |
---|
[258] | 544 | parameter. This subsection works in exactly the same way as the pixel |
---|
[819] | 545 | subsection discussed in \S\ref{sec-input}. The \texttt{StatSec} will |
---|
| 546 | be trimmed if necessary so that it lies wholly within the image |
---|
| 547 | subsection being used (\ie that given by the \texttt{subsection} |
---|
| 548 | parameter - this governs what pixels are read in and so are able to be |
---|
| 549 | used in the calculations). |
---|
[258] | 550 | |
---|
[819] | 551 | Note that \texttt{StatSec} applies only to the statistics used to |
---|
| 552 | determine the threshold. It does not affect the calculation of |
---|
| 553 | statistics in the case of the wavelet reconstruction. Note also that |
---|
| 554 | pixels flagged as BLANK or as part of the ``Milky Way'' range of |
---|
| 555 | channels are ignored in the statistics calculations. |
---|
| 556 | |
---|
[258] | 557 | \secC{Determining the threshold} |
---|
| 558 | |
---|
| 559 | Once the statistics have been calculated, the threshold is determined |
---|
| 560 | in one of two ways. The first way is a simple sigma-clipping, where a |
---|
| 561 | threshold is set at a fixed number $n$ of standard deviations above |
---|
| 562 | the mean, and pixels above this threshold are flagged as detected. The |
---|
[386] | 563 | value of $n$ is set with the parameter \texttt{snrCut}. The ``mean'' |
---|
| 564 | and ``standard deviation'' here are estimated according to |
---|
| 565 | \texttt{flagRobustStats}, as discussed in \S\ref{sec-stats}. In this |
---|
| 566 | first case only, if the user specifies a threshold, using the |
---|
| 567 | \texttt{threshold} parameter, the sigma-clipped value is ignored. |
---|
[258] | 568 | |
---|
[158] | 569 | The second method uses the False Discovery Rate (FDR) technique |
---|
| 570 | \citep{miller01,hopkins02}, whose basis we briefly detail here. The |
---|
| 571 | false discovery rate (given by the number of false detections divided |
---|
| 572 | by the total number of detections) is fixed at a certain value |
---|
| 573 | $\alpha$ (\eg $\alpha=0.05$ implies 5\% of detections are false |
---|
| 574 | positives). In practice, an $\alpha$ value is chosen, and the ensemble |
---|
| 575 | average FDR (\ie $\langle FDR \rangle$) when the method is used will |
---|
| 576 | be less than $\alpha$. One calculates $p$ -- the probability, |
---|
| 577 | assuming the null hypothesis is true, of obtaining a test statistic as |
---|
| 578 | extreme as the pixel value (the observed test statistic) -- for each |
---|
| 579 | pixel, and sorts them in increasing order. One then calculates $d$ |
---|
| 580 | where |
---|
| 581 | \[ |
---|
| 582 | d = \max_j \left\{ j : P_j < \frac{j\alpha}{c_N N} \right\}, |
---|
| 583 | \] |
---|
| 584 | and then rejects all hypotheses whose $p$-values are less than or |
---|
| 585 | equal to $P_d$. (So a $P_i<P_d$ will be rejected even if $P_i \geq |
---|
| 586 | j\alpha/c_N N$.) Note that ``reject hypothesis'' here means ``accept |
---|
| 587 | the pixel as an object pixel'' (\ie we are rejecting the null |
---|
| 588 | hypothesis that the pixel belongs to the background). |
---|
| 589 | |
---|
[277] | 590 | The $c_N$ value here is a normalisation constant that depends on the |
---|
[158] | 591 | correlated nature of the pixel values. If all the pixels are |
---|
| 592 | uncorrelated, then $c_N=1$. If $N$ pixels are correlated, then their |
---|
| 593 | tests will be dependent on each other, and so $c_N = \sum_{i=1}^N |
---|
| 594 | i^{-1}$. \citet{hopkins02} consider real radio data, where the pixels |
---|
[265] | 595 | are correlated over the beam. For the calculations done in \duchamp, |
---|
[801] | 596 | $N = B \times C$, where $B$ is the beam area in pixels, calculated |
---|
[571] | 597 | from the FITS header (if the correct keywords -- BMAJ, BMIN -- are not |
---|
[801] | 598 | present, the size of the beam is taken from the input parameters - see |
---|
| 599 | discussion in \S\ref{sec-results}, and if these parameters are not |
---|
| 600 | given, $B=1$), and $C$ is the number of neighbouring channels that can |
---|
| 601 | be considered to be correlated. |
---|
[158] | 602 | |
---|
[543] | 603 | The use of the FDR method is governed by the \texttt{flagFDR} flag, |
---|
| 604 | which is \texttt{false} by default. To set the relevant parameters, |
---|
| 605 | use \texttt{alphaFDR} to set the $\alpha$ value, and |
---|
| 606 | \texttt{FDRnumCorChan} to set the $C$ value discussed above. These |
---|
| 607 | have default values of 0.01 and 2 respectively. |
---|
| 608 | |
---|
[158] | 609 | The theory behind the FDR method implies a direct connection between |
---|
| 610 | the choice of $\alpha$ and the fraction of detections that will be |
---|
[265] | 611 | false positives. These detections, however, are individual pixels, |
---|
| 612 | which undergo a process of merging and rejection (\S\ref{sec-merger}), |
---|
| 613 | and so the fraction of the final list of detected objects that are |
---|
| 614 | false positives will be much smaller than $\alpha$. See the discussion |
---|
| 615 | in \S\ref{sec-notes}. |
---|
[158] | 616 | |
---|
[265] | 617 | %\secC{Storage of detected objects in memory} |
---|
| 618 | % |
---|
| 619 | %It is useful to understand how \duchamp stores the detected objects in |
---|
| 620 | %memory while it is running. This makes use of nested C++ classes, so |
---|
| 621 | %that an object is stored as a class that includes the set of detected |
---|
| 622 | %pixels, plus all the various calculated parameters (fluxes, WCS |
---|
| 623 | %coordinates, pixel centres and extrema, flags,...). The set of pixels |
---|
| 624 | %are stored using another class, that stores 3-dimensional objects as a |
---|
| 625 | %set of channel maps, each consisting of a $z$-value and a |
---|
| 626 | %2-dimensional object (a spatial map if you like). This 2-dimensional |
---|
| 627 | %object is recorded using ``run-length'' encoding, where each row (a |
---|
| 628 | %fixed $y$ value) is stored by the starting $x$-value and the length |
---|
[158] | 629 | |
---|
[691] | 630 | \secB{Merging, growing and rejecting detected objects} |
---|
[158] | 631 | \label{sec-merger} |
---|
| 632 | |
---|
[691] | 633 | \secC{Merging} |
---|
[158] | 634 | |
---|
[819] | 635 | The searches described above are either 1- or 2-dimensional only. They |
---|
[691] | 636 | do not know anything about the third dimension that is likely to be |
---|
[819] | 637 | present. To build up 3D sources, merging of detections must be |
---|
[691] | 638 | done. This is done via an algorithm that matches objects judged to be |
---|
| 639 | ``close'', according to one of two criteria. |
---|
| 640 | |
---|
[158] | 641 | One criterion is to define two thresholds -- one spatial and one in |
---|
| 642 | velocity -- and say that two objects should be merged if there is at |
---|
| 643 | least one pair of pixels that lie within these threshold distances of |
---|
| 644 | each other. These thresholds are specified by the parameters |
---|
| 645 | \texttt{threshSpatial} and \texttt{threshVelocity} (in units of pixels |
---|
| 646 | and channels respectively). |
---|
| 647 | |
---|
| 648 | Alternatively, the spatial requirement can be changed to say that |
---|
| 649 | there must be a pair of pixels that are \emph{adjacent} -- a stricter, |
---|
| 650 | but perhaps more realistic requirement, particularly when the spatial |
---|
[691] | 651 | pixels have a large angular size (as is the case for \hi |
---|
| 652 | surveys). This method can be selected by setting the parameter |
---|
| 653 | \texttt{flagAdjacent=true} in the parameter file. The velocity |
---|
| 654 | thresholding is always done with the \texttt{threshVelocity} test. |
---|
[158] | 655 | |
---|
[691] | 656 | |
---|
| 657 | \secC{Stages of merging} |
---|
| 658 | |
---|
| 659 | This merging can be done in two stages. The default behaviour is for |
---|
| 660 | each new detection to be compared with those sources already detected, |
---|
| 661 | and for it to be merged with the first one judged to be close. No |
---|
| 662 | other examination of the list is done at this point. |
---|
| 663 | |
---|
| 664 | This step can be turned off by setting |
---|
| 665 | \texttt{flagTwoStageMerging=false}, so that new detections are simply |
---|
| 666 | added to the end of the list, leaving all merging to be done in the |
---|
| 667 | second stage. |
---|
| 668 | |
---|
| 669 | The second, main stage of merging is more thorough, Once the searching |
---|
| 670 | is completed, the list is iterated through, looking at each pair of |
---|
| 671 | objects, and merging appropriately. The merged objects are then |
---|
| 672 | included in the examination, to see if a merged pair is suitably close |
---|
| 673 | to a third. |
---|
| 674 | |
---|
| 675 | \secC{Growing} |
---|
| 676 | |
---|
[158] | 677 | Once the detections have been merged, they may be ``grown''. This is a |
---|
[265] | 678 | process of increasing the size of the detection by adding nearby |
---|
| 679 | pixels (according to the \texttt{threshSpatial} and |
---|
| 680 | \texttt{threshVelocity} parameters) that are above some secondary |
---|
[766] | 681 | threshold and not already part of a detected object. This threshold |
---|
| 682 | should be lower than the one used for the initial detection, but above |
---|
| 683 | the noise level, so that faint pixels are only detected when they are |
---|
| 684 | close to a bright pixel. This threshold is specified via one of two |
---|
| 685 | input parameters. It can be given in terms of the noise statistics via |
---|
| 686 | \texttt{growthCut} (which has a default value of $3\sigma$), or it can |
---|
| 687 | be directly given via \texttt{growthThreshold}. Note that if you have |
---|
| 688 | given the detection threshold with the \texttt{threshold} parameter, |
---|
| 689 | the growth threshold \textbf{must} be given with |
---|
| 690 | \texttt{growthThreshold}. |
---|
[265] | 691 | |
---|
| 692 | The use of the growth algorithm is controlled by the |
---|
[158] | 693 | \texttt{flagGrowth} parameter -- the default value of which is |
---|
| 694 | \texttt{false}. If the detections are grown, they are sent through the |
---|
[766] | 695 | merging algorithm a second time, to pick up any detections that should |
---|
| 696 | be merged at the new lower threshold (\ie they have grown into each |
---|
| 697 | other). |
---|
[158] | 698 | |
---|
[691] | 699 | \secC{Rejecting} |
---|
| 700 | \label{sec-reject} |
---|
| 701 | |
---|
[158] | 702 | Finally, to be accepted, the detections must span \emph{both} a |
---|
[285] | 703 | minimum number of channels (enabling the removal of any spurious |
---|
| 704 | single-channel spikes that may be present), and a minimum number of |
---|
| 705 | spatial pixels. These numbers, as for the original detection step, are |
---|
[720] | 706 | set with the \texttt{minChannels}, \texttt{minPix} parameters and |
---|
| 707 | \texttt{minVoxels}. The channel requirement means a source must have |
---|
| 708 | at least one set of \texttt{minChannels} consecutive channels to be |
---|
| 709 | accepted. The spatial pixels (\texttt{minPix}) requirement refers to |
---|
| 710 | distinct spatial pixels (which are possibly in different channels), |
---|
| 711 | while the voxels requirement refers to the total number of voxels |
---|
| 712 | detected. |
---|
[482] | 713 | |
---|
[691] | 714 | It is possible to do this rejection stage before the main merging and |
---|
| 715 | growing stage. This could be done to remove narrow (hopefully |
---|
| 716 | spurious) sources from the list before growing them, to reduce the |
---|
| 717 | number of false positives in the final list. This mode can be selected |
---|
[696] | 718 | by setting the input parameter \texttt{flagRejectBeforeMerge=true} - |
---|
| 719 | caution is urged if you use this in conjunction with |
---|
| 720 | \texttt{flagTwoStageMerging=false}. |
---|