Source code for lentil.util

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

[docs]def circle(shape, radius, shift=(0, 0)): """Compute a circle with anti-aliasing. Parameters ---------- shape : array_like Size of output in pixels (nrows, ncols) radius : float Radius of circle in pixels 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[0]) + np.square(cc - shift[1])) return np.clip(radius + 0.5 - r, 0.0, 1.0)
[docs]def circlemask(shape, radius, shift=(0, 0)): """Compute a circular mask. Parameters ---------- shape : array_like Size of output in pixels (nrows, ncols) radius : float Radius of circle in pixels shift : array_like, optional How far to shift center in float (rows, cols). Default is (0, 0). Returns ------- mask : ndarray """ mask =, radius, shift) mask[mask > 0] = 1 return mask
[docs]def hexagon(shape, radius, rotate=False): """Compute a hexagon mask. Parameters ---------- shape : array_like Size of output in pixels (nrows, ncols) radius : int 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 ------- mask : ndarray """ inner_radius = radius * np.sqrt(3)/2 side_length = radius/2 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))) mask = np.zeros(shape) mask[rect] = 1 mask[left_tri] = 1 mask[right_tri] = 1 if rotate: return mask.transpose() else: return mask
[docs]def slit(shape, width, length=None): """Compute a slit mask. 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 ------- mask : ndarray """ rr, cc = lentil.helper.mesh(shape) slit = np.ones(shape) length = shape[1] 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 =, img.ravel()) c =, img.ravel()) return r, c
[docs]def pad(array, shape): """Zero-pad an array. Note that pad works for both two and three dimensional arrays. Parameters ---------- array : array_like Array to be padded. shape : tuple of ints Shape of output array in ``(nrows, ncols)``. Returns ------- padded : ndarray 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[0] - array.shape[0+offset] dc = shape[1] - array.shape[1+offset] if dr <= 0: rmin0 = (array.shape[0+offset] - shape[0])//2 rmax0 = rmin0 + shape[0] rmin1 = 0 rmax1 = shape[0] else: rmin0 = 0 rmax0 = array.shape[0+offset] rmin1 = (shape[0] - array.shape[0+offset])//2 rmax1 = rmin1 + array.shape[0+offset] if dc <= 0: cmin0 = (array.shape[1+offset] - shape[1])//2 cmax0 = cmin0 + shape[1] cmin1 = 0 cmax1 = shape[1] else: cmin0 = 0 cmax0 = array.shape[1] cmin1 = (shape[1] - array.shape[1+offset])//2 cmax1 = cmin1 + array.shape[1+offset] if array.ndim < 3: padded = np.zeros((shape[0], shape[1]), dtype=array.dtype) padded[rmin1:rmax1, cmin1:cmax1] = array[rmin0:rmax0, cmin0:cmax0] else: padded = np.zeros((array.shape[0], shape[0], shape[1]), dtype=array.dtype) padded[:, rmin1:rmax1, cmin1:cmax1] = array[:, rmin0:rmax0, cmin0:cmax0] return padded
[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), 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[0] = (slice[1] - slice[0]) = (r_end - r_start) shape[1] = (slice[3] - slice[2]) = (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[1] - slice[0]) == shape[0] assert(slice[3] - slice[2]) == shape[1] # return the requested view. Note that numpy will implicitly handle a # third dimension if one is present return img[slice[0]:slice[1], slice[2]:slice[3]] else: # return the padded array. Note that pad will handle a third dimension # if one exists return lentil.pad(img, shape)
[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][[0, -1]] cmin, cmax = np.where(cols)[0][[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 See Also -------- :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[0], img.shape[1]//factor, img.shape[2]//factor) img_rebinned = np.zeros(rebinned_shape, dtype=img.dtype) for i in range(img.shape[0]): img_rebinned[i] = img[i].reshape(rebinned_shape[1], factor, rebinned_shape[2], factor).sum(-1).sum(1) else: img_rebinned = img.reshape(img.shape[0]//factor, factor, img.shape[1]//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. mask : array_like, optional Binary mask applied after rescaling. If None (default), a mask is 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. See Also -------- :func:`rebin` """ img = np.asarray(img) if mask is None: # take the real portion to ensure that even if img is complex, mask will # be real mask = np.zeros_like(img).real mask[img != 0] = 1 if shape is None: shape = np.ceil((img.shape[0]*scale, img.shape[1]*scale)).astype(int) else: if np.isscalar(shape): shape = np.ceil((shape*scale, shape*scale)).astype(int) else: shape = np.ceil((shape[0]*scale, shape[1]*scale)).astype(int) x = (np.arange(shape[1], dtype=np.float64) - shape[1]/2.)/scale + img.shape[1]/2. y = (np.arange(shape[0], dtype=np.float64) - shape[0]/2.)/scale + img.shape[0]/2. xx, yy = np.meshgrid(x, y) mask = map_coordinates(mask, [yy, xx], order=1, mode='nearest') mask[mask < np.finfo(mask.dtype).eps] = 0 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) out *= mask 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[0] * npix[0]), num/(min_q * du[1] * npix[1])
[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