1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
#!/usr/bin/env python3
#
# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT License
#
"""
Taper the sky image (input of OSKAR simulation) to mitigate the
side lobes effects, which causes trouble in creating good images.
The circular Tukey window is adopted, which is also built in by
e.g., WSClean.
"""
import os
import sys
import argparse
import logging
from datetime import datetime, timezone
import numpy as np
from scipy import interpolate
import astropy.io.fits as fits
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(os.path.basename(sys.argv[0]))
def make2d(w1d, x=None):
"""
Create 2D filter from the 1D one.
Parameters
----------
w1d : 1D `~numpy.ndarray`
The input 1D filter/window
x : 1D `~numpy.ndarray`, optional
The X-axis values of the input ``w1d`` filter
Returns
-------
w2d : 2D `~numpy.ndarray`
Created 2D filter/window from the input 1D one.
Credit
------
[1] MATLAB - creating 2D convolution filters
https://cn.mathworks.com/matlabcentral/newsreader/view_thread/23588
"""
if x is None:
L = len(w1d)
M = (L-1) / 2
x = np.linspace(-M, M, num=L)
xmax = np.max(x)
xsize = int(2 * xmax + 1)
xx = np.linspace(-xmax, xmax, num=xsize)
xg, yg = np.meshgrid(xx, xx)
r = np.sqrt(xg**2 + yg**2)
ridx = (r <= xmax)
w2d = np.zeros(shape=(xsize, xsize))
finterp = interpolate.interp1d(x=x, y=w1d, kind="linear")
w2d[ridx] = finterp(r[ridx])
return w2d
def tukey(L, Minner=0.5, Mouter=0.5, Lmin=0):
"""
The 1D Tukey filter/window
^
| Lmin | Minner*L | constant ones | Mouter*L |
+--------------------------------------------. (L)
Parameters
----------
L : int
The length (number of points) of the window/filter
Minner : float
The fraction of ``L`` for the inner (left), increasing part
Mouter : float
The fraction of ``L`` for the outer (left), decreasing part
Lmin : int
The beginning point for the inner (left) part
References
----------
* MATLAB - Tukey (tapered cosine) window
https://cn.mathworks.com/help/signal/ref/tukeywin.html
* WSClean - Weight tapering
https://sourceforge.net/p/wsclean/wiki/Tapering/
"""
w = np.zeros(L)
Linner = int(Minner * L)
Louter = int(Mouter * L)
Lones = L - Lmin - Linner - Louter
if Lones < 0:
raise RuntimeError("invalid parameters combination")
w[(Lmin+Linner):(Lmin+Linner+Lones)] = 1.0
if Linner > 0:
xinner = np.arange(0, Linner)
x2 = (np.pi/Linner) * (xinner - Linner)
w[Lmin:(Lmin+Linner)] = 0.5 + 0.5 * np.cos(x2)
if Louter > 0:
xouter = np.arange(0, Louter)
x2 = (np.pi/Louter) * xouter
w[(Lmin+Linner+Lones):] = 0.5 + 0.5 * np.cos(x2)
return w
def main():
parser = argparse.ArgumentParser(
description="Taper sky image with Tukey window")
parser.add_argument("-C", "--clobber", dest="clobber",
action="store_true",
help="overwrite existing file")
parser.add_argument("-r", "--window-inner", dest="Rinner",
type=int, required=True,
help="inner radius [pixel] of the tukey window")
parser.add_argument("-R", "--window-outer", dest="Router",
type=int, required=True,
help="outer radius [pixel] of the tukey window")
parser.add_argument("-T", "--outfile-taper", dest="outfile_taper",
help="save the applied taper into a FITS image")
parser.add_argument("infile", help="input FITS sky image")
parser.add_argument("outfile", help="output tapered sky image")
args = parser.parse_args()
if os.path.exists(args.outfile):
if args.clobber:
logger.warning("Removed existing file: %s" % args.outfile)
os.remove(args.outfile)
else:
raise OSError("Output file already exists: %s" % args.outfile)
if args.outfile_taper and os.path.exists(args.outfile_taper):
if args.clobber:
logger.warning("Removed existing file: %s" % args.outfile_taper)
os.remove(args.outfile_taper)
else:
raise OSError("Output file already exists: %s" %
args.outfile_taper)
with fits.open(args.infile) as f:
image = f[0].data.astype(np.float32)
header = f[0].header.copy(strip=True)
logger.info("Read sky image from file: %s" % args.infile)
L = args.Router
Louter = L - args.Rinner
Mouter = Louter / L
logger.info("Create Tukey window: Rin=%d, Rout=%d" %
(args.Rinner, args.Router))
w1d = tukey(L, Minner=0, Mouter=Mouter)
x1d = np.arange(len(w1d))
logger.info("Make 2D Tukey window")
w2d = make2d(w1d, x1d)
xw = w2d.shape[0]
ximg = image.shape[0]
if xw == ximg:
taper = w2d
elif xw > ximg:
i1 = int((xw-ximg) / 2)
i2 = int((xw+ximg) / 2)
taper = w2d[i1:i2, i1:i2]
else:
taper = np.zeros(shape=image.shape)
i1 = int((ximg-xw) / 2)
i2 = int((ximg+xw) / 2)
taper[i1:i2, i1:i2] = w2d
logger.info("Taper the sky image ...")
image *= taper
header["DATE"] = (datetime.now(timezone.utc).astimezone().isoformat(),
"File creation date")
header.add_history(" ".join(sys.argv))
hdu = fits.PrimaryHDU(data=image, header=header)
hdu.writeto(args.outfile)
logger.info("Wrote tapered sky image to file: %s" % args.outfile)
if args.outfile_taper:
header["OBJECT"] = ("Tukey Window", "Taper window")
header["Rinner"] = (args.Rinner, "[pixel] inner radius")
header["Router"] = (args.Router, "[pixel] outer radius")
hdu = fits.PrimaryHDU(data=taper.astype(np.float32), header=header)
hdu.writeto(args.outfile_taper)
logger.info("Wrote tapered to file: %s" % args.outfile_taper)
if __name__ == "__main__":
main()
|
