| 1 | import numpy
|
|---|
| 2 |
|
|---|
| 3 | def image(himage, psf, dirty, niter, thresh):
|
|---|
| 4 | himage.open(psf)
|
|---|
| 5 | psfdata = himage.getchunk()
|
|---|
| 6 | print 'psf read, max val %s' % numpy.max(psfdata)
|
|---|
| 7 | psfdata /= numpy.max(psfdata)
|
|---|
| 8 | print 'psf normalised'
|
|---|
| 9 |
|
|---|
| 10 | # if not dirty:
|
|---|
| 11 | # sky = numpy.zeros(numpy.shape(psf))
|
|---|
| 12 | # dirty = numpy.zeros(numpy.shape(psf))
|
|---|
| 13 |
|
|---|
| 14 | # sky[70,70]=5000
|
|---|
| 15 | # for i in range(numpy.shape(psf)[0]):
|
|---|
| 16 | # for j in range(numpy.shape(psf)[1]):
|
|---|
| 17 | # if pow((i-20),2) + pow((j-20),2) < 100:
|
|---|
| 18 | # sky[i,j] = 10000
|
|---|
| 19 |
|
|---|
| 20 | # for i in range(numpy.shape(psf)[0]):
|
|---|
| 21 | # for j in range(numpy.shape(psf)[1]):
|
|---|
| 22 | # a1o, a2o=overlapIndices(sky, psf,
|
|---|
| 23 | # i+1-sky.shape[0]/2,
|
|---|
| 24 | # j+1-sky.shape[1]/2)
|
|---|
| 25 | # dirty[a1o[0]:a1o[1],a1o[2]:a1o[3]]+=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*sky[i,j]
|
|---|
| 26 | # else:
|
|---|
| 27 | # himage.open(dirty)
|
|---|
| 28 | # dirtydata = himage.getchunk()
|
|---|
| 29 | # himage.fromarray(outfile='dirty.im', pixels=dirty, overwrite=True)
|
|---|
| 30 | # print 'dirty read'
|
|---|
| 31 |
|
|---|
| 32 | himage.open(dirty)
|
|---|
| 33 | dirtydata = himage.getchunk()
|
|---|
| 34 | print 'dirty read'
|
|---|
| 35 |
|
|---|
| 36 | comps, res = hogbom(dirtydata, psfdata, True, 0.05, thresh, niter)
|
|---|
| 37 | print 'clean completed'
|
|---|
| 38 |
|
|---|
| 39 | himage.fromarray(outfile='residual.im', pixels=res, overwrite=True)
|
|---|
| 40 | himage.fromarray(outfile='comps.im', pixels=comps, overwrite=True)
|
|---|
| 41 | himage.fromarray(outfile='image.im', pixels=comps+res, overwrite=True)
|
|---|
| 42 | himage.close()
|
|---|
| 43 |
|
|---|
| 44 |
|
|---|
| 45 | # Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
|
|---|
| 46 | # Initial version August 2010
|
|---|
| 47 | #
|
|---|
| 48 | # This file is part of pydeconv. This work is licensed under GNU GPL
|
|---|
| 49 | # V2 (http://www.gnu.org/licenses/gpl.html)
|
|---|
| 50 | """
|
|---|
| 51 | Clean based deconvolution, using numpy
|
|---|
| 52 | """
|
|---|
| 53 |
|
|---|
| 54 | def overlapIndices(a1, a2,
|
|---|
| 55 | shiftx, shifty):
|
|---|
| 56 | if shiftx >=0:
|
|---|
| 57 | a1xbeg=shiftx
|
|---|
| 58 | a2xbeg=0
|
|---|
| 59 | a1xend=a1.shape[0]
|
|---|
| 60 | a2xend=a1.shape[0]-shiftx
|
|---|
| 61 | else:
|
|---|
| 62 | a1xbeg=0
|
|---|
| 63 | a2xbeg=-shiftx
|
|---|
| 64 | a1xend=a1.shape[0]+shiftx
|
|---|
| 65 | a2xend=a1.shape[0]
|
|---|
| 66 |
|
|---|
| 67 | if shifty >=0:
|
|---|
| 68 | a1ybeg=shifty
|
|---|
| 69 | a2ybeg=0
|
|---|
| 70 | a1yend=a1.shape[1]
|
|---|
| 71 | a2yend=a1.shape[1]-shifty
|
|---|
| 72 | else:
|
|---|
| 73 | a1ybeg=0
|
|---|
| 74 | a2ybeg=-shifty
|
|---|
| 75 | a1yend=a1.shape[1]+shifty
|
|---|
| 76 | a2yend=a1.shape[1]
|
|---|
| 77 |
|
|---|
| 78 | return (a1xbeg, a1xend, a1ybeg, a1yend), (a2xbeg, a2xend, a2ybeg, a2yend)
|
|---|
| 79 |
|
|---|
| 80 |
|
|---|
| 81 |
|
|---|
| 82 | def hogbom(dirty,
|
|---|
| 83 | psf,
|
|---|
| 84 | window,
|
|---|
| 85 | gain,
|
|---|
| 86 | thresh,
|
|---|
| 87 | niter):
|
|---|
| 88 | """
|
|---|
| 89 | Hogbom clean
|
|---|
| 90 |
|
|---|
| 91 | :param dirty: The dirty image, i.e., the image to be deconvolved
|
|---|
| 92 |
|
|---|
| 93 | :param psf: The point spread-function
|
|---|
| 94 |
|
|---|
| 95 | :param window: Regions where clean components are allowed. If
|
|---|
| 96 | True, thank all of the dirty image is assumed to be allowed for
|
|---|
| 97 | clean components
|
|---|
| 98 |
|
|---|
| 99 | :param gain: The "loop gain", i.e., the fraction of the brightest
|
|---|
| 100 | pixel that is removed in each iteration
|
|---|
| 101 |
|
|---|
| 102 | :param thresh: Cleaning stops when the maximum of the absolute
|
|---|
| 103 | deviation of the residual is less than this value
|
|---|
| 104 |
|
|---|
| 105 | :param niter: Maximum number of components to make if the
|
|---|
| 106 | threshold "thresh" is not hit
|
|---|
| 107 | """
|
|---|
| 108 | comps=numpy.zeros(dirty.shape)
|
|---|
| 109 | res=numpy.array(dirty)
|
|---|
| 110 | if window is True:
|
|---|
| 111 | window=numpy.ones(dirty.shape,
|
|---|
| 112 | numpy.bool)
|
|---|
| 113 | for i in range(niter):
|
|---|
| 114 | mx, my=numpy.unravel_index(numpy.fabs(res[window]).argmax(), res.shape)
|
|---|
| 115 | print 'mx, my', mx, my
|
|---|
| 116 | mval=res[mx, my]*gain
|
|---|
| 117 | print 'mval', mval
|
|---|
| 118 | comps[mx, my]+=mval
|
|---|
| 119 | a1o, a2o=overlapIndices(dirty, psf,
|
|---|
| 120 | mx+1-dirty.shape[0]/2,
|
|---|
| 121 | my+1-dirty.shape[1]/2)
|
|---|
| 122 | # mx-dirty.shape[0]/2,
|
|---|
| 123 | # my-dirty.shape[1]/2)
|
|---|
| 124 | print 'a1o', a1o
|
|---|
| 125 | print 'a2o', a2o
|
|---|
| 126 | res[a1o[0]:a1o[1],a1o[2]:a1o[3]]-=psf[a2o[0]:a2o[1],a2o[2]:a2o[3]]*mval
|
|---|
| 127 | print '59,56', res[59,56]
|
|---|
| 128 | if numpy.fabs(res).max() < thresh:
|
|---|
| 129 | break
|
|---|
| 130 | return comps, res
|
|---|
