filter_fchop.py
3.16 KB
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
title = "Fchop"
tip = "chops high frequencies to fft data"
onein = True
import numpy as np
import scipy.signal as sps
from guidata.qt.QtGui import QMessageBox
from guidata.dataset.datatypes import DataSet
from guidata.dataset.dataitems import (IntItem, FloatArrayItem, StringItem,
ChoiceItem, FloatItem, DictItem,
BoolItem)
from guiqwt.config import _
def exp2(x):
if x > -12.0: y = np.exp(x)
else: y = 0.0
return y
class NOD3_App:
def __init__(self, parent):
self.parent = parent
self.parent.activateWindow()
def Error(self, msg):
QMessageBox.critical(self.parent.parent(), title,
_(u"Error:")+"\n%s" % str(msg))
def compute_app(self):
class Param(DataSet):
rchop = FloatItem('RChop', default=0.5, min=0.0, max=1.0)
#xbeam = FloatItem('XBeam', default=0.0, min=0.0, max=100.0)
#ybeam = FloatItem('YBeam', default=0.0, min=0.0, max=100.0)
#param = Param(_("Plait"), "Set parameters")
param = None
name = title.replace(" ", "")
self.parent.compute_11(name, lambda m, p: self.function(m, p), param, onein)
def ftgrd2(self, data, header, rchop, xbeam2, ybeam2):
#rows, cols = data.shape
cols, rows = data.shape
#xsamp, ysamp = abs(header['CDELT1']), abs(header['CDELT2'])
ysamp, xsamp = abs(header['CDELT1']), abs(header['CDELT2'])
if 'BMAJ' in header:
xbeam1 = header['BMAJ']
else:
xbeam1 = 2.5*xsamp
if 'BMIN' in header:
ybeam1 = header['BMIN']
else:
ybeam1= 2.5*ysamp
if 'DIAMETER' in header:
d = header['DIAMETER']
else:
d = 100.0
if 'CTYPE3' in header:
if header['CTYPE3'] == 'LAMBDA':
rlam = header['CRVAL3']
else:
rlam = 299.792458/float(header['CRVAL3'])*1000000.0
else:
rlam = 0.1
if not 'BMAJ' in header:
self.Error("BMAJ / BMIN not in header")
return None
#rchop = xsamp/xbeam1
spfreq = 2./min(xbeam1, ybeam1) / 1.22
fchopsq = abs(rchop)*spfreq**2
f1xsq = (1.0/(xsamp * float(cols-1)))**2
f1ysq = (1.0/(ysamp * float(rows-1)))**2
x = np.arange(int(cols))
y = np.arange(int(rows))
xx = x - ((x+0)/int(cols/2+2))*cols
yy = y - ((y+0)/int(rows/2+2))*rows
rxsq = np.outer(np.ones(cols), yy*yy*f1ysq)
rysq = np.outer(xx*xx*f1xsq, np.ones(rows))
data = np.where(rxsq+rysq > fchopsq, 0.0 + 0j, data)
return data
def function(self, m, p):
mask = np.isnan(m.data)
xmin = np.nanmin(m.data)
m.data = np.where(np.isnan(m.data), xmin, m.data)
fftdata = np.fft.fft2(m.data)
#fftdata = self.ftgrd2(fftdata, m.header, p.rchop, 0.0, 0.0)
fftdata = self.ftgrd2(fftdata, m.header, 0.333333, 0.0, 0.0)
if type(fftdata) == type(None):
return [], p
m.data = np.fft.ifft2(fftdata).real
m.data = np.where(mask, np.nan, m.data)
return m, p