pynkowski.data.array

  1import numpy as np
  2from .base_da import DataField
  3
  4
  5
  6def _eval_field(field, pix):
  7    return field[tuple(pix[:,ii] for ii in np.arange(field.ndim))]
  8
  9def _hotspot_array(field):
 10    """Find the local maxima of the input field.
 11
 12    Arguments
 13    ---------
 14    field : np.array
 15        Array of the pixelized field values.
 16
 17    Returns
 18    -------
 19    pixels : np.array
 20        Indices of the pixels which are local maxima.
 21
 22    values : np.array
 23        Values of input map which are local maxima.
 24
 25    """
 26    ndim = len(field.shape)
 27    # First we shift the field in each direction by 1 pixel, and check that the original pixel is larger than all the `2**ndim` neighbours.
 28    max_mask = np.all(field > np.array([np.roll(field, shift, axis=ii) for shift in [-1,1] for ii in range(ndim)]), axis=0)
 29
 30    # We then remove all the pixels in the border to remove the edge effects.
 31    for dim in range(max_mask.ndim):
 32        # Make current dimension the first dimension
 33        array_moved = np.moveaxis(max_mask, dim, 0)
 34        # Set border values to `False` in the current dimension
 35        array_moved[0] = False
 36        array_moved[-1] = False
 37        # No need to reorder the dimensions as moveaxis returns a view
 38    max_mask = np.pad(max_mask, pad_width=1, mode='constant', constant_values=False)
 39
 40    pixels = np.argwhere(max_mask)
 41    values = _eval_field(field, pixels)
 42    return(pixels, values)
 43
 44class DataArray(DataField):
 45    """Class for a pixelized Euclidean data cube.
 46    
 47    Parameters
 48    ----------
 49    field : np.array
 50        Pixelized field as a 3D array. All pixels (or voxels) are expected to be the same size. The pixel shape is not necessarily a cube (see `spacing`).
 51        The field can have different sizes in each dimension.
 52        
 53    normalise : bool, optional
 54        If `True`, the map is normalised to have unit variance. Note: the mean of the map is not forced to be 0.
 55        Default: `True`.
 56        
 57    mask : np.array or None, optional
 58        Mask where the field if considered. It is a bool array of the same shape that `field`.
 59        Default: all data is included.
 60
 61    spacing : float or np.array, optional
 62        Spacing between pixels (centres) in each dimension. If a float, the spacing is the same in all dimensions. 
 63        If an array, it must have the same length as the number of dimensions of the field.
 64        
 65    Attributes
 66    ----------
 67    field : np.array
 68        Data of the field as a 3D array.
 69        
 70    dim : int
 71        Dimension of the space where the field is defined. In this case, the space is the 3D space and this is 3.
 72        
 73    name : str
 74        Name of the field. In this case, "data cube"
 75        
 76    first_der : np.array or None
 77        First **covariant** derivatives of the field in an orthonormal basis of the space. Same structure as `field`, and shape `(dim, field.shape)`.
 78        
 79    second_der : np.array or None
 80        Second **covariant** derivatives of the field in an orthonormal basis of the space. Same structure as `field`, and shape `(dim*(dim+1)/2, field.shape)`.
 81
 82    spacing : float or np.array
 83        Spacing between pixels (centres) in each dimension. If a float, the spacing is the same in all dimensions.
 84        If an array, it must have the same length as the number of dimensions of the field.
 85        
 86    """
 87    def __init__(self, field, normalise=True, mask=None, spacing=1.):
 88        dim = len(field.shape)
 89        super().__init__(field, dim=dim, name='DataArray', mask=mask)
 90        if self.mask.shape != self.field.shape:
 91            raise ValueError('The map and the mask have different shapes')
 92        if normalise:
 93            self.field /= np.sqrt(self.get_variance())
 94        self.first_der = None
 95        self.second_der = None
 96        self.spacing = spacing
 97            
 98    def get_variance(self):
 99        """Compute the variance of the field within the mask. 
100
101        Returns
102        -------
103        var : float
104            The variance of the input field within the mask.
105
106        """    
107        return (np.var(self.field[self.mask]))
108
109    def get_first_der(self):
110        """Compute the covariant first derivatives of the field. 
111        It stores the derivatives in order. E.g., in a 3D array:
112        
113        - first covariant derivative wrt e₁ in self.first_der[0]
114        - first covariant derivative wrt e₂ in self.first_der[1]
115        - first covariant derivative wrt e₃ in self.first_der[2]
116
117        """
118        self.first_der = np.array(np.gradient(self.field, self.spacing, edge_order=2))
119
120    def get_second_der(self):
121        """Compute the covariant second derivatives of the field. 
122        It stores the second derivatives in the following order: first, all the diagonal derivatives, then the cross derivatives with e₁, then with e₂, and so on.
123        E.g., in a 3D array:
124        
125        - second covariant derivative wrt e₁e₁ in self.second_der[0]
126        - second covariant derivative wrt e₂e₂ in self.second_der[1]
127        - second covariant derivative wrt e₃e₃ in self.second_der[2]
128        - second covariant derivative wrt e₁e₂ in self.second_der[3]
129        - second covariant derivative wrt e₁e₃ in self.second_der[4]
130        - second covariant derivative wrt e₂e₃ in self.second_der[5]
131
132        """
133        if self.first_der is None:
134            self.get_first_der()
135        self.second_der = np.zeros((self.dim*(self.dim+1)//2, *self.field.shape))
136        d = self.dim
137        for i in np.arange(self.dim, dtype=int):
138            indeces = [i]
139            for j in np.arange(i+1, self.dim, dtype=int):
140                indeces.append((d*(d-1)/2) - (d-i)*((d-i)-1)/2 + j - i - 1 + d)
141            indeces = np.array(indeces, dtype=int)
142            self.second_der[indeces] = np.gradient(self.first_der[i], self.spacing, edge_order=2)[i:d]
143
144    def maxima_list(self):
145        """Find the local maxima of the field.
146
147        Returns
148        -------
149        values : np.array
150            Values of input map which are local maxima.
151
152        """
153        pixels, values = _hotspot_array(self.field)
154        return values[_eval_field(self.mask, pixels)]
155    
156    def minima_list(self):
157        """Find the local minima of the field.
158
159        Returns
160        -------
161        values : np.array
162            Values of input map which are local minima.
163
164        """
165        pixels, values = _hotspot_array(-self.field)
166        return -values[_eval_field(self.mask, pixels)]
167
168
169
170__all__ = ["DataArray"]
171__docformat__ = "numpy"