SpheMakh / HI-Inator

Radio interefometry simulator/imager tailored for HI sky models.
GNU General Public License v2.0
1 stars 0 forks source link

Add (iterative) cleaning using masks #6

Open SpheMakh opened 9 years ago

SpheMakh commented 9 years ago

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:

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 commented 9 years ago

@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 commented 9 years ago

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

SpheMakh commented 9 years ago

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 commented 9 years ago

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