# Source code for lentil.util

import numpy as np
from scipy.ndimage import map_coordinates
import lentil

"""Compute a circle with anti-aliasing.

Parameters
----------
shape : array_like
Size of output in pixels (nrows, ncols)

shift : (2,) array_like, optional
How far to shift center in float (rows, cols). Default is (0, 0).

Returns
-------
circle : ndarray

"""
rr, cc = lentil.helper.mesh(shape)
r = np.sqrt(np.square(rr - shift) + np.square(cc - shift))
return np.clip(radius + 0.5 - r, 0.0, 1.0)

Parameters
----------
shape : array_like
Size of output in pixels (nrows, ncols)

shift : array_like, optional
How far to shift center in float (rows, cols). Default is (0, 0).

Returns
-------

"""

Parameters
----------
shape : array_like
Size of output in pixels (nrows, ncols)

Radius of outscribing circle (which also equals the side length) in
pixels.

rotate : bool
Rotate mask so that flat sides are aligned with the Y direction instead
of the default orientation which is aligned with the X direction.

Returns
-------

"""

rr, cc = lentil.helper.mesh(shape)

rect = np.where((np.abs(cc) <= side_length) & (np.abs(rr) <= inner_radius))
left_tri = np.where((cc <= -side_length) & (cc >= -radius) & (np.abs(rr) <= (cc + radius)*np.sqrt(3)))
right_tri = np.where((cc >= side_length) & (cc <= radius) & (np.abs(rr) <= (radius - cc)*np.sqrt(3)))

if rotate:
else:

[docs]def slit(shape, width, length=None):

Parameters
----------
shape : array_like
Size of output in pixels (nrows, ncols)
width : float
Slit width in pixels
length : float, optional
Slit length in pixels. If not specified, the slit spans
the entire column shape (default).

Returns
-------

"""

rr, cc = lentil.helper.mesh(shape)
slit = np.ones(shape)

length = shape if length is None else length
width_clip = np.clip(0.5 + (width/2) - np.abs(rr), 0, 1)
length_clip = np.clip(0.5 + (length/2) - np.abs(cc), 0, 1)

slit = np.minimum(np.minimum(slit, width_clip), length_clip)

return slit

[docs]def centroid(img):
"""Compute image centroid location.

Parameters
----------
img : array_like
Input array.

Returns
-------
tuple
(r,c) centroid location.
"""

img = np.asarray(img)
img = img/np.sum(img)
nr, nc = img.shape
rr, cc = np.mgrid[0:nr, 0:nc]

r = np.dot(rr.ravel(), img.ravel())
c = np.dot(cc.ravel(), img.ravel())

return r, c

Note that pad works for both two and three dimensional arrays.

Parameters
----------
array : array_like

shape : tuple of ints
Shape of output array in (nrows, ncols).

Returns
-------
Zero-padded array with shape (nrows, ncols). If array has a
third dimension, the return shape will be (nrows, ncols, depth).
"""

array = np.asarray(array)

offset = 0
if array.ndim == 3:
offset = 1

dr = shape - array.shape[0+offset]
dc = shape - array.shape[1+offset]

if dr <= 0:
rmin0 = (array.shape[0+offset] - shape)//2
rmax0 = rmin0 + shape
rmin1 = 0
rmax1 = shape
else:
rmin0 = 0
rmax0 = array.shape[0+offset]
rmin1 = (shape - array.shape[0+offset])//2
rmax1 = rmin1 + array.shape[0+offset]

if dc <= 0:
cmin0 = (array.shape[1+offset] - shape)//2
cmax0 = cmin0 + shape
cmin1 = 0
cmax1 = shape
else:
cmin0 = 0
cmax0 = array.shape
cmin1 = (shape - array.shape[1+offset])//2
cmax1 = cmin1 + array.shape[1+offset]

if array.ndim < 3:
else:
padded = np.zeros((array.shape, shape, shape), dtype=array.dtype)
padded[:, rmin1:rmax1, cmin1:cmax1] = array[:, rmin0:rmax0, cmin0:cmax0]

[docs]def window(img, shape=None, slice=None):
"""Extract an appropriately sized, potentially windowed array

Parameters
----------
img : array_like
Data to window

shape : array_like or None, optional
Output shape given as (nrows, ncols). If None (default), an
uncropped array is returned.

slice : array_like or None, optional
Indices of img array to return given as (r_start, r_end, c_start,
c_end). This definition follows standard numpy indexing.

Returns
-------
data : ndarray
Trimmed and windowed img array

Notes
-----
* If img is a single value (img.size == 1), self.data is returned
regardless of what shape and slice are.
* If shape is given but slice is None, the returned ndarray
is trimmed about the center of the array using :func:lentil.pad.
* If slice is given but shape is None, the returned ndarray
is extracted from img according to the indices in slice
* If both shape and slice are given, the returned ndarray is
extracted from img according to the indices in slice and the
following expressions must also be true:

.. code::

shape = (slice - slice) = (r_end - r_start)
shape = (slice - slice) = (c_end - c_start)

"""

img = np.asarray(img)
if img.size == 1:
return img

if shape is None and slice is None:
return img
elif slice is not None:
if shape is not None:
# ensure size consistency
assert(slice - slice) == shape
assert(slice - slice) == shape

# return the requested view. Note that numpy will implicitly handle a
# third dimension if one is present
return img[slice:slice, slice:slice]

else:
# return the padded array. Note that pad will handle a third dimension
# if one exists

[docs]def boundary(x, threshold=0):
"""Find bounding row and column indices of data within an array.

Parameters
----------
x : array_like
Input array

threshold : float, optional
Masking threshold to apply before boundary finding. Only values
in x that are larger than threshold are considered in the boundary
finding operation. Default is 0.

Returns
-------
rmin, rmax, cmin, cmax : ints
Boundary indices

"""
x = np.asarray(x)
x = (x > threshold)

rows = np.any(x, axis=1)
cols = np.any(x, axis=0)

rmin, rmax = np.where(rows)[[0, -1]]
cmin, cmax = np.where(cols)[[0, -1]]

return rmin, rmax, cmin, cmax

[docs]def rebin(img, factor):
"""Rebin an image by an integer factor.

Parameters
----------
img : array_like
Array or cube of arrays to rebin. If a cube is provided, the first dimension
should index the image slices.

factor : int
Rebinning factor

Returns
-------
img : ndarray
Rebinned image

--------
:func:rescale

"""
img = np.asarray(img)

if np.iscomplexobj(img):
raise ValueError('rebin is not defined for complex data')

if img.ndim == 3:
rebinned_shape = (img.shape, img.shape//factor, img.shape//factor)
img_rebinned = np.zeros(rebinned_shape, dtype=img.dtype)
for i in range(img.shape):
img_rebinned[i] = img[i].reshape(rebinned_shape, factor,
rebinned_shape, factor).sum(-1).sum(1)
else:
img_rebinned = img.reshape(img.shape//factor, factor, img.shape//factor,
factor).sum(-1).sum(1)

return img_rebinned

[docs]def rescale(img, scale, shape=None, mask=None, order=3, mode='nearest',
unitary=True):
"""Rescale an image by interpolation.

Parameters
----------
img : array_like
Image to rescale

scale : float
Scaling factor. Scale factors less than 1 will shrink the image. Scale
factors greater than 1 will grow the image.

shape : array_like or int, optional
Output shape. If None (default), the output shape will be the input img
shape multiplied by the scale factor.

created from the nonzero portions of img. To skip masking operation,
set mask = np.ones_like(img)

order : int, optional
Order of spline interpolation used for rescaling operation. Default is
3. Order must be in the range 0-5.

mode : {'constant', 'nearest', 'reflect', 'wrap'}, optional
Points outside the boundaries of the input are filled according to the
given mode. Default is 'constant'.

unitary : bool, optional
Normalization flag. If True (default), a normalization is performed on
the output such that the rescaling operation is unitary and image power
(if complex) or intensity (if real) is conserved.

Returns
-------
img : ndarray
Rescaled image.

Note
----
The post-rescale masking operation should have no real effect on the
resulting image but is included to eliminate interpolation artifacts that
sometimes appear in large clusters of zeros in rescaled images.

--------
:func:rebin

"""

img = np.asarray(img)

# take the real portion to ensure that even if img is complex, mask will
# be real

if shape is None:
shape = np.ceil((img.shape*scale, img.shape*scale)).astype(int)
else:
if np.isscalar(shape):
shape = np.ceil((shape*scale, shape*scale)).astype(int)
else:
shape = np.ceil((shape*scale, shape*scale)).astype(int)

x = (np.arange(shape, dtype=np.float64) - shape/2.)/scale + img.shape/2.
y = (np.arange(shape, dtype=np.float64) - shape/2.)/scale + img.shape/2.

xx, yy = np.meshgrid(x, y)

if np.iscomplexobj(img):
out = np.zeros(shape, dtype=np.complex128)
out.real = map_coordinates(img.real, [yy, xx], order=order, mode=mode)
out.imag = map_coordinates(img.imag, [yy, xx], order=order, mode=mode)
else:
out = map_coordinates(img, [yy, xx], order=order, mode=mode)

if unitary:
out *= np.sum(img)/np.sum(out)

return out

[docs]def pixelscale_nyquist(wave, f_number):
"""Compute the output plane sampling which is Nyquist sampled for
intensity.

Parameters
----------
wave : float
Wavelength in meters

f_number : float
Optical system F/#

Returns
-------
float
Sampling in meters

"""
return f_number * wave / 2

[docs]def min_sampling(wave, z, du, npix, min_q):
num = np.min(wave) * z
return num/(min_q * du * npix), num/(min_q * du * npix)

[docs]def normalize_power(array, power=1):
r"""Normalizie the power in an array.

The total power in an array is given by

.. math::

P = \sum{\left|\mbox{array}\right|^2}

A normalization coefficient is computed as

.. math::

c = \sqrt{\frac{p}{\sum{\left|\mbox{array}\right|^2}}}

The array returned by a will be scaled by the normalization coefficient so
that its power is equal to :math:p.

Parameters
----------
array : array_like
Array to be normalized

power : float, optional
Desired power in normalized array. Default is 1.

Returns
-------
array : ndarray
Normalized array

"""
array = np.asarray(array)
return array * np.sqrt(power/np.sum(np.abs(array)**2))

[docs]def sanitize_shape(shape):
shape = np.asarray(shape)
if shape.size == 0:
shape = ()
else:
if shape.shape == ():
shape = np.append(shape, shape)
return tuple(shape)

[docs]def sanitize_bandpass(vec):
vec = np.asarray(vec)
if vec.shape == ():
vec = vec[np.newaxis, ...]
return vec