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