Field Operations¶

def extend(
    self,
    x=(0, 0),
    y=(0, 0),
    z=(0, 0),
    boundary_x=("periodic", "periodic"),
    boundary_y=("periodic", "periodic"),
    boundary_z=("periodic", "periodic"),
) -> "Field":
    """Return a field extended using MICRESS-style boundary conditions.

    Args:
        x (tuple of int, optional): Lower and upper extension in the x
            direction. Defaults to ``(0, 0)``.
        y (tuple of int, optional): Lower and upper extension in the y
            direction. Defaults to ``(0, 0)``.
        z (tuple of int, optional): Lower and upper extension in the z
            direction. Defaults to ``(0, 0)``.
        boundary_x (tuple of str, optional): Lower and upper boundary
            conditions for x. Defaults to ``("periodic", "periodic")``.
        boundary_y (tuple of str, optional): Lower and upper boundary
            conditions for y. Defaults to ``("periodic", "periodic")``.
        boundary_z (tuple of str, optional): Lower and upper boundary
            conditions for z. Defaults to ``("periodic", "periodic")``.

            ``"insulating"`` mirrors while duplicating the outermost
            sample, ``"symmetric"`` mirrors without duplicating it, and
            ``"periodic"`` wraps. Extensions and boundary pairs are
            specified in physical x/y/z order.

    Returns:
        Field: Extended field with the original time and spacing.

    Raises:
        ValueError: If the field is not three-dimensional, an extension is
            negative, a boundary condition is invalid, or periodic
            conditions are not used on both sides of an axis.
    """
    if self.ndim != 3:
        raise ValueError("field must be 3D (nz, ny, nx)")

    data = self._extend_array(
        np.asarray(self),
        extension_zyx=(z, y, x),
        boundary_zyx=(boundary_z, boundary_y, boundary_x),
    )
    return Field(data, time=self.time, spacing=self.spacing)

def resample(
    self,
    shape: Tuple[int, int, int],
    method: str = "nearest",
) -> "Field":
    """Resample the field to a new cell geometry.

    The physical extent of the field is preserved. Consequently, changing
    the number of cells changes the spacing accordingly.

    Args:
        shape (Tuple[int, int, int]): Target shape in NumPy axis order,
            i.e. ``(z, y, x)``.
        method (str, optional): Resampling method. Defaults to ``"nearest"``.
            ``"nearest"`` samples the nearest source cell, preserving the
            input dtype and suiting categorical data such as grain or phase
            IDs. ``"linear"`` interpolates between source cell centres and
            returns floating-point data for continuous fields. ``"mean"``
            computes an overlap-weighted mean of source cells covered by
            each target cell, returning floating-point data for conservative
            coarsening of continuous cell data.

    Returns:
        Field: Resampled field with updated spacing and preserved time.

    Raises:
        ValueError: If the target shape or method is invalid, or if
            ``"mean"`` would refine an axis.
    """
    if self.ndim != 3:
        raise ValueError("resample() requires a three-dimensional Field")

    try:
        shape = tuple(int(n) for n in shape)
    except (TypeError, ValueError) as exc:
        raise ValueError("shape must contain three positive integers") from exc

    if len(shape) != 3 or any(n <= 0 for n in shape):
        raise ValueError("shape must contain three positive integers")

    method = method.lower()
    if method not in {"nearest", "linear", "mean"}:
        raise ValueError("method must be 'nearest', 'linear', or 'mean'")

    old_shape = tuple(self.shape)

    if method == "mean" and any(new > old for old, new in zip(old_shape, shape)):
        raise ValueError(
            "'mean' can only be used for coarsening; "
            "use 'nearest' or 'linear' when refining a field"
        )

    data = np.asarray(self)

    if shape == old_shape:
        result = data.copy()

    elif method == "nearest":
        result = data

        for axis, (old_n, new_n) in enumerate(zip(old_shape, shape)):
            if old_n == new_n:
                continue

            # Map target cell centres into source-cell coordinates.
            coords = (
                np.arange(new_n, dtype=np.float64) + 0.5
            ) * old_n / new_n - 0.5

            indices = np.floor(coords + 0.5).astype(np.intp)
            indices = np.clip(indices, 0, old_n - 1)

            result = np.take(result, indices, axis=axis)

    elif method == "linear":
        result = data.astype(np.result_type(data.dtype, np.float64), copy=False)

        for axis, (old_n, new_n) in enumerate(zip(old_shape, shape)):
            if old_n == new_n:
                continue

            # A singleton axis contains no information from which to
            # interpolate. Repeating the sole value is the natural result.
            if old_n == 1:
                result = np.repeat(result, new_n, axis=axis)
                continue

            # Coordinates of target cell centres expressed in source
            # cell-centre index coordinates.
            coords = (
                np.arange(new_n, dtype=np.float64) + 0.5
            ) * old_n / new_n - 0.5

            # Constant extrapolation at the physical boundaries.
            coords = np.clip(coords, 0.0, old_n - 1.0)

            lower = np.floor(coords).astype(np.intp)
            upper = np.minimum(lower + 1, old_n - 1)
            weight = coords - lower

            lower_values = np.take(result, lower, axis=axis)
            upper_values = np.take(result, upper, axis=axis)

            weight_shape = [1] * result.ndim
            weight_shape[axis] = new_n
            weight = weight.reshape(weight_shape)

            result = lower_values * (1.0 - weight) + upper_values * weight

    else:  # method == "mean"
        result = data.astype(np.result_type(data.dtype, np.float64), copy=False)

        for axis, (old_n, new_n) in enumerate(zip(old_shape, shape)):
            if old_n == new_n:
                continue

            # Source cells occupy intervals [i, i + 1]. Target cells cover
            # the same physical interval [0, old_n], divided into new_n
            # equal cells.
            target_edges = np.linspace(0.0, float(old_n), new_n + 1)

            source_left = np.arange(old_n, dtype=np.float64)
            source_right = source_left + 1.0

            target_left = target_edges[:-1, None]
            target_right = target_edges[1:, None]

            overlap = np.maximum(
                0.0,
                np.minimum(target_right, source_right)
                - np.maximum(target_left, source_left),
            )

            # Divide by target-cell width so every row sums to one.
            weights = overlap / (float(old_n) / new_n)

            # Move the current spatial axis to the front, apply the
            # one-dimensional overlap operator, then restore the axis.
            moved = np.moveaxis(result, axis, 0)
            moved_shape = moved.shape

            result = weights @ moved.reshape(old_n, -1)
            result = result.reshape((new_n,) + moved_shape[1:])
            result = np.moveaxis(result, 0, axis)

    spacing = None
    if self.spacing is not None:
        spacing = tuple(
            old_spacing * old_n / new_n
            for old_spacing, old_n, new_n in zip(self.spacing, old_shape, shape)
        )

    return Field(result, time=self.time, spacing=spacing)