Skip to content

Commit

Permalink
Ostrich (#15)
Browse files Browse the repository at this point in the history
* Ignoring invalid division by zero or nan
An ostrich approach to #11

* Dropzweights (#14)

* remove full nan slices

* don't acr by default

* look for max of model only at unflagged locations

* drop invalid gausspars

* accommodate different keys for different header

* idx -> fidx

* fresq -> freqs

* Exclude any bands that might be awful with `-ds`

* pep8 indetation level

* syntax error fix

* Update spi_fitter.py

* Update spi_fitter.py

---------

Co-authored-by: landmanbester <[email protected]>

* Fix freq dimension check
Fix empty gauspar checking

---------

Co-authored-by: Athanaseus Javas Ramaila <[email protected]>
Co-authored-by: landmanbester <[email protected]>
  • Loading branch information
3 people authored Apr 10, 2024
1 parent a22c45d commit fef25e9
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 10 deletions.
19 changes: 12 additions & 7 deletions spimple/apps/image_convolver.py
Original file line number Diff line number Diff line change
Expand Up @@ -63,14 +63,18 @@ def image_convolver():
l_coord -= ref_l
m_coord, ref_m = data_from_header(hdr, axis=2)
m_coord -= ref_m
if hdr["CTYPE4"].lower() == 'freq':
freq_axis = 4
elif hdr["CTYPE3"].lower() == 'freq':
freq_axis = 3

if hdr["NAXIS"] >2:
if "CTYPE4" in hdr and hdr["CTYPE4"].lower() == 'freq':
freq_axis = 4
elif hdr["CTYPE3"].lower() == 'freq':
freq_axis = 3
else:
raise ValueError("Freq axis must be 3rd or 4th")
freqs, ref_freq = data_from_header(hdr, axis=freq_axis)
else:
raise ValueError("Freq axis must be 3rd or 4th")
freqs, ref_freq = data_from_header(hdr, axis=freq_axis)

freqs = np.array([0])

nchan = freqs.size
gausspari = ()
if freqs.size > 1:
Expand Down Expand Up @@ -151,6 +155,7 @@ def image_convolver():
if imagei.ndim != 3:
raise ValueError("Unsupported number of image dimensions")
print('Convolving image', file=log)

image, gausskern = convolve2gaussres(imagei, xx, yy, gaussparf,
opts.nthreads, gausspari,
opts.padding_frac)
Expand Down
18 changes: 15 additions & 3 deletions spimple/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -89,8 +89,16 @@ def Gaussian2D(xin, yin, GaussPar=(1., 1., 0.), normalise=True):
S0, S1, PA = GaussPar
Smaj = np.maximum(S0, S1)
Smin = np.minimum(S0, S1)
A = np.array([[1. / Smin ** 2, 0],

# hackaround for runtime warning
if 0 in (Smaj, Smin):
print("Nan values have been encountered. Subverting RuntimeWarning")
A = np.array([[np.nan, 0],
[0, np.nan]])
else:
A = np.array([[1. / Smin ** 2, 0],
[0, 1. / Smaj ** 2]])


c, s, t = np.cos, np.sin, np.deg2rad(-PA)
R = np.array([[c(t), -s(t)],
Expand Down Expand Up @@ -157,15 +165,19 @@ def convolve2gaussres(image, xx, yy, gaussparf, nthreads, gausspari=None, pfrac=
imhat = r2c(iFs(image, axes=ax), axes=ax, forward=True, nthreads=nthreads, inorm=0)

# convolve to desired resolution
if gausspari is None:
if gausspari in [None, ()]:
imhat *= gausskernhat
else:
for i in range(nband):
thiskern = Gaussian2D(xx, yy, gausspari[i], normalise=norm_kernel).astype(np.float64)
thiskern = np.pad(thiskern[None], padding, mode='constant')
thiskernhat = r2c(iFs(thiskern, axes=ax), axes=ax, forward=True, nthreads=nthreads, inorm=0)

convkernhat = np.where(np.abs(thiskernhat)>1e-10, gausskernhat/thiskernhat, 0.0)
if not np.all(np.isnan(thiskernhat)):
convkernhat = np.where(np.abs(thiskernhat)>1e-10, gausskernhat/thiskernhat, 0.0)
else:
print("Nan values have been encountered. Subverting RuntimeWarning")
convkernhat = np.zeros_like(thiskernhat).astype("complex")

imhat[i] *= convkernhat[0]

Expand Down

0 comments on commit fef25e9

Please sign in to comment.