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
|
---|