| 1 | #import numpy
|
|---|
| 2 |
|
|---|
| 3 | # Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
|
|---|
| 4 | # Initial version August 2010
|
|---|
| 5 | #
|
|---|
| 6 | # This file is part of pydeconv. This work is licensed under GNU GPL
|
|---|
| 7 | # V2 (http://www.gnu.org/licenses/gpl.html)
|
|---|
| 8 | """
|
|---|
| 9 | Clean based deconvolution, using numpy
|
|---|
| 10 | """
|
|---|
| 11 |
|
|---|
| 12 | # Note added by jfl. I think a1 and a2 are assumed the same shape.
|
|---|
| 13 |
|
|---|
| 14 | def overlapIndices(a1, a2,
|
|---|
| 15 | shiftx, shifty):
|
|---|
| 16 | if shiftx >=0:
|
|---|
| 17 | a1xbeg=shiftx
|
|---|
| 18 | a2xbeg=0
|
|---|
| 19 | a1xend=a1.shape[0]
|
|---|
| 20 | a2xend=a1.shape[0]-shiftx
|
|---|
| 21 | else:
|
|---|
| 22 | a1xbeg=0
|
|---|
| 23 | a2xbeg=-shiftx
|
|---|
| 24 | a1xend=a1.shape[0]+shiftx
|
|---|
| 25 | a2xend=a1.shape[0]
|
|---|
| 26 |
|
|---|
| 27 | if shifty >=0:
|
|---|
| 28 | a1ybeg=shifty
|
|---|
| 29 | a2ybeg=0
|
|---|
| 30 | a1yend=a1.shape[1]
|
|---|
| 31 | a2yend=a1.shape[1]-shifty
|
|---|
| 32 | else:
|
|---|
| 33 | a1ybeg=0
|
|---|
| 34 | a2ybeg=-shifty
|
|---|
| 35 | a1yend=a1.shape[1]+shifty
|
|---|
| 36 | a2yend=a1.shape[1]
|
|---|
| 37 |
|
|---|
| 38 | return (a1xbeg, a1xend, a1ybeg, a1yend), (a2xbeg, a2xend, a2ybeg, a2yend)
|
|---|
| 39 |
|
|---|
| 40 |
|
|---|
| 41 |
|
|---|
| 42 | def hogbom(dirty,
|
|---|
| 43 | psf,
|
|---|
| 44 | window,
|
|---|
| 45 | gain,
|
|---|
| 46 | thresh,
|
|---|
| 47 | niter):
|
|---|
| 48 | """
|
|---|
| 49 | Hogbom clean
|
|---|
| 50 |
|
|---|
| 51 | :param dirty: The dirty image, i.e., the image to be deconvolved
|
|---|
| 52 |
|
|---|
| 53 | :param psf: The point spread-function
|
|---|
| 54 |
|
|---|
| 55 | :param window: Regions where clean components are allowed. If
|
|---|
| 56 | True, thank all of the dirty image is assumed to be allowed for
|
|---|
| 57 | clean components
|
|---|
| 58 |
|
|---|
| 59 | :param gain: The "loop gain", i.e., the fraction of the brightest
|
|---|
| 60 | pixel that is removed in each iteration
|
|---|
| 61 |
|
|---|
| 62 | :param thresh: Cleaning stops when the maximum of the absolute
|
|---|
| 63 | deviation of the residual is less than this value
|
|---|
| 64 |
|
|---|
| 65 | :param niter: Maximum number of components to make if the
|
|---|
| 66 | threshold "thresh" is not hit
|
|---|
| 67 | """
|
|---|
| 68 | comps=numpy.zeros(dirty.shape)
|
|---|
| 69 | res=numpy.array(dirty)
|
|---|
| 70 | grid=numpy.indices(dirty.shape)
|
|---|
| 71 | if window is True:
|
|---|
| 72 | window = [slice(dirty.shape), slice(dirty.shape)]
|
|---|
| 73 |
|
|---|
| 74 | for i in range(niter):
|
|---|
| 75 | argmax=numpy.unravel_index(numpy.fabs(res[window]).argmax(), res[window].shape)
|
|---|
| 76 | mx=grid[0][window][argmax]
|
|---|
| 77 | my=grid[1][window][argmax]
|
|---|
| 78 |
|
|---|
| 79 | mval=res[mx, my]*gain
|
|---|
| 80 | comps[mx, my]+=mval
|
|---|
| 81 | a1o, a2o=pythonclean.overlapIndices(dirty, psf,
|
|---|
| 82 | mx-dirty.shape[0]/2,
|
|---|
| 83 | my-dirty.shape[1]/2)
|
|---|
| 84 | res[a1o[0]:a1o[1],a1o[2]:a1o[3]]-=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*mval
|
|---|
| 85 | if numpy.fabs(res).max() < thresh:
|
|---|
| 86 | break
|
|---|
| 87 | return comps, res
|
|---|
| 88 |
|
|---|
