Add (iterative) cleaning using masks

[SpheMakh]

This is how it is approximately done in the miriad pipeline used for
Serra et al. 2012 http://esoads.eso.org/abs/2012MNRAS.422.1835S or
Wang et al. 2013 http://esoads.eso.org/abs/2013MNRAS.433..270W or
Janowiecki et al. 2015 http://esoads.eso.org/abs/2015ApJ...801...96J
which is HI work. HI can be concentrated, point source-like, or extended. The trick is to filter the cubes with different size filters (3d Gaussians) to enhance emission at different scales and to define the clean region based on that.

Use a Hogbom or Clarke clean (certainly other variants work as well). Number of iterations infinite, but with a cutoff level set.

Repeat the loop below niters times with the running index of the iteration being denoted as i . (0 <= i < niters)

For the ith iteration define a number cleancut[i] with cleancut[i] < cleancut[i+i] for all i and the last cleancut[niters-1] being very large.

For the ith iteration define a multiplicator n[i] for the rms noise to determine clean masks, n[i]>=n[i+1].

For the ith iteration define a cutoff paramter cutoff[i] for the cleaning, cutoff[i] >= cutoff[i+1].

Start:

• Define clean region:
  -- produce very few similar cubes by filtering the restored cube from the previous iteration (dirty cube at the beginning) with m different 3D Gaussians (with FWHMs fwhm_x_1,fwhm_y_1, fwhm_v_1, ..., fwhm_x_m, fwhm_y_m, fwhm_v_m), in all three directions. (Remark: in the pipeline we performed Hanning smoothing along the velocity axis instead, just because Miriad does not provide Gaussian filtering on the third axis.)
  -- For each resulting cube and for the original cube calculate the maximum max and the rms.
  -- For each resulting cube and for the original cube calculate a threshold thre as the maximum of max/cleancut[i] and n[i]*rms.
  -- For each resulting cube and for the original cube define a partial clean region as all voxels above thre
  -- Define the united clean region as the union of the clean regions of the single cubes
• Clean and restore:
  -- Clean the dirty cube using the united clean region, i.e. allowing only for clean components inside the united clean region. Use an infinite possible number of iterations, but set a cutoff parameter of cutoff[i]*rms_original, where rms_original is the rms of the original (unfiltered) data cube
  -- Produce the restored cube.
• Check if loop is repeated:
  -- if i = niters stop, otherwise increase i by 1 and start the loop again, goto Start

For the Serra pipeline the parameters are as follows:

niters = 4
cleancut = [3, 6, 9, 1000]
n = [5, 5, 5, 5]
cutoff = [1, 1, 1, 1]
m = 1
fwhm_x_1 = fwhm_y_1 = 60''
fwhm_v_1 ~ 5 Pixels

(In that pipeline a Hanning filter with a width of 9 pixels was used rather than a Gaussian.) A code snippet might clarify some quesitons:

cleancut=[3,6,9,1000]
for jj in range(len(cleancut)):    
    CONVOL('r'+str(jj),60,'r'+str(jj)+'_60',han=9)
    Max=IMAGESTAT('r'+str(jj),domax=1)
    Sig=IMAGESTAT('r'+str(jj),dosig=1)
    Max60=IMAGESTAT('r'+str(jj)+'_60,domax=1)
    Sig60=IMAGESTAT('r'+str(jj)+'_60,dosig=1)
    Thr=str(max(float(Max)/cleancut[jj],5*float(Sig)))
    Thr60=str(max(float(Max60)/cleancut[jj],5*float(Sig60)))
    MAKEMASK('r0,'r'+str(jj)+'.gt.'+Thr+'.or.r'+str(jj)+'_60'+'.gt.'+Thr60)
    CLEAN('r0','b,'c'+str(jj),ctf=1*Sig)

Remarks:
i) To determine the rms (and filtering out all positive sources), I prefer the following method: Make a histogram and fit a Gaussian to the left part from the peak (left indicating the negative direction), neglecting (ideally) the positive part of the histogram. The sigma of that Gaussian is a good estimate of the rms.
ii) The Serra pipeline runs only one filter to determine an additional clean mask of a convolved cube. For WSRT data it might be good to also try a

fwhm_x_1 = fwhm_y_1 = 30''
fwhm_v_1 ~ 3 Pixels

in addition.
iii) Again, this does not need to be restricted to a Hogbom clean, but can be used with any suitable deconvolution algorithm.

[SpheMakh]

@gigjozsa, from this statement: "cutoff[i]*rms_original, where rms_original is the rms of the original (unfiltered) data cube"

Is the "original unfiltered cube" the dirty map, or is it a clean cube without without filtering?

[gigjozsa]

Yup, that's too cryptic. It's the restored data cube in this iteration, not convolved (and in that sense the "original").

[SpheMakh]

I see. Do we use the noise per channel or for the whole cube? I suppose since we are cleaning each channel it makes sense to use the noise per channel?

[gigjozsa]

Both not perfect. But per channel appears to be more appropriate. Any highly variable noise would be a concern anyway.