Module curlew.geometry

Utility functions for generating basic geometries (grids, sections, etc.) and performing other simple geometric tasks.

Functions

def blended_wave(x, f=0.5, A=1, T=6.283185307179586)

Expand source code

def blended_wave( x, f=0.5, A=1, T=2*np.pi ):
    """
    Blend a triangle and a sinusoid wave to get varying degrees of sharpness in the fold geometry. 

    Parameters:
    x : array-like
        Input values (time or position).
    f : float
        The blending factor to use. 0 gives a perfectly sinusoidal wave, while 1 gives a perfectly triangular wave.
    A : float
        Amplitude of the wave.
    T : float
        Period (wavelength) of the function.

    Returns:
    f_x : array-like
        Approximated blend between a triangle wave and a sine wave. 
    """
    y1 = triangle_wave(x, A=A, T=T, n_terms=1) # sinusoid wave 
    y2 = triangle_wave(x, A=A, T=T, n_terms=11) # triangle wave 
    return (1-f)*y1 + f*y2

Blend a triangle and a sinusoid wave to get varying degrees of sharpness in the fold geometry.

Parameters: x : array-like Input values (time or position). f : float The blending factor to use. 0 gives a perfectly sinusoidal wave, while 1 gives a perfectly triangular wave. A : float Amplitude of the wave. T : float Period (wavelength) of the function.

Returns: f_x : array-like Approximated blend between a triangle wave and a sine wave.

def compute_topology_adjacency(soft_masks, symmetrize=True, connection='conservative')

Expand source code

def compute_topology_adjacency(soft_masks, symmetrize=True, connection="conservative"):
    """
    Differentiable adjacency matrix from class soft-assignment volumes. This is
    used by the `curlew.core.Geode.topology` method to compute the adjacency matrix from
    gridded lithology and structure IDs.

    Uses depthwise convolution with a neighbourhood kernel, then batch matrix
    multiplication to accumulate voxel contact area between classes.

    Parameters
    ----------
    soft_masks : torch.Tensor
        Shape ``(B, C, *spatial)`` with ``spatial`` either ``(H, W)`` or ``(D, H, W)``.
        Values are non-negative class weights (e.g. softmax probabilities) per voxel.
    symmetrize : bool, optional
        If True (default), return ``(A + A.T) / 2`` so contact area is mutual.
    connection : str, optional
        ``'conservative'`` (default): only face-shared neighbours count (4-connected in 2D,
        6-connected in 3D). ``'inclusive'``: face, edge, and vertex neighbours count
        (8-connected in 2D, 26-connected in 3D).

    Returns
    -------
    torch.Tensor
        Shape ``(B, C, C)``; diagonal entries are zero (no self-contact).
    """
    # Require batched multi-channel image/volume layout: (B, C, H, W) or (B, C, D, H, W).
    if soft_masks.ndim not in (4, 5):
        raise ValueError(
            "soft_masks must have shape (B, C, H, W) or (B, C, D, H, W); "
            f"got {tuple(soft_masks.shape)}"
        )
    # Unpack batch size, class count, and voxel count for the matrix multiply below.
    B, C = soft_masks.shape[:2]
    spatial = soft_masks.shape[2:]
    n_vox = int(np.prod(spatial))
    device, dtype = soft_masks.device, soft_masks.dtype

    # Conservative = shared faces only; inclusive = shared faces, edges, or vertices.
    connection = connection.lower()
    if connection not in ("conservative", "inclusive"):
        raise ValueError(
            f"connection must be 'conservative' or 'inclusive', not {connection!r}"
        )

    # Build a depthwise stencil: centre 0, ones where two cells are considered adjacent.
    if len(spatial) == 2:
        if connection == "conservative":
            # Face neighbours only (4-connected).
            kernel = torch.tensor(
                [[0.0, 1.0, 0.0], [1.0, 0.0, 1.0], [0.0, 1.0, 0.0]],
                device=device,
                dtype=dtype,
            )
        else:
            # Face, edge, and vertex neighbours (8-connected).
            kernel = torch.ones((3, 3), device=device, dtype=dtype)
            kernel[1, 1] = 0.0
        kernel = kernel.view(1, 1, 3, 3).expand(C, 1, 3, 3)
        dilated = F.conv2d(soft_masks, kernel, padding=1, groups=C)
    elif len(spatial) == 3:
        kernel = torch.zeros((3, 3, 3), device=device, dtype=dtype)
        if connection == "conservative":
            # Face neighbours only (6-connected).
            kernel[1, 1, 0] = kernel[1, 1, 2] = 1.0
            kernel[1, 0, 1] = kernel[1, 2, 1] = 1.0
            kernel[0, 1, 1] = kernel[2, 1, 1] = 1.0
        else:
            # Face, edge, and vertex neighbours (26-connected).
            kernel[:] = 1.0
            kernel[1, 1, 1] = 0.0
        kernel = kernel.view(1, 1, 3, 3, 3).expand(C, 1, 3, 3, 3)
        dilated = F.conv3d(soft_masks, kernel, padding=1, groups=C)
    else:
        raise ValueError(f"topology supports 2D and 3D grids only; spatial shape {spatial}")

    # Collapse spatial dims so bmm accumulates contact area: A[i,j] = sum_v X_i(v) * Neigh(X_j)(v).
    flat_masks = soft_masks.reshape(B, C, n_vox)
    flat_dilated = dilated.reshape(B, C, n_vox).transpose(1, 2)
    adjacency = torch.bmm(flat_masks, flat_dilated)
    
    # Remove self-contact (every voxel neighbours its own class).
    adjacency = adjacency * (1.0 - torch.eye(C, device=device, dtype=dtype).unsqueeze(0))
    
    # Average with transpose so edge voxels do not bias i->j vs j->i.
    if symmetrize:
        adjacency = 0.5 * (adjacency + adjacency.transpose(1, 2))
    return adjacency

Differentiable adjacency matrix from class soft-assignment volumes. This is used by the Geode.topology() method to compute the adjacency matrix from gridded lithology and structure IDs.

Uses depthwise convolution with a neighbourhood kernel, then batch matrix multiplication to accumulate voxel contact area between classes.

Parameters

soft_masks : torch.Tensor

Shape (B, C, *spatial) with spatial either (H, W) or (D, H, W). Values are non-negative class weights (e.g. softmax probabilities) per voxel.

symmetrize : bool, optional

If True (default), return (A + A.T) / 2 so contact area is mutual.

connection : str, optional

'conservative' (default): only face-shared neighbours count (4-connected in 2D, 6-connected in 3D). 'inclusive': face, edge, and vertex neighbours count (8-connected in 2D, 26-connected in 3D).

Returns

torch.Tensor

Shape (B, C, C); diagonal entries are zero (no self-contact).

def grid(extent: tuple | numpy.ndarray, size: tuple, **kwds)

Expand source code

def grid( extent : tuple | np.ndarray, size : tuple, **kwds ):
    """
    Construct a grid using the provided extent (or the bounds of passed point set) and
    voxel size (float) or grid dimensions (int). 
    
    Parameters
    ----------
    extent : tuple | np.ndarray
        The extent of the grid in each dimension as [(xmin,ymin,...), (xmax, ymax, ...)]. If a tuple is provided, it is interpreted as the bounds of the grid in each dimension.
        If a numpy array is provided, the min and max of these points in each dimension will be computed and used as the bounds.
    size : tuple
        The size of the grid in each dimension. If a float is provided, it is interpreted as the voxel size in each dimension.
        If an integer is provided, it is interpreted as the number of grid cells in each dimension.
    **kwds : dict
        Additional keyword arguments are passed directly to `Grid.__init__`.
    """
    # get grid bounds if not passed directly
    if isinstance(extent, np.ndarray):
        extent = (extent.min(axis=0), extent.max(axis=0))

    ndim = len(extent[0])
    if isinstance(size, (float, int)):
        size = tuple([size] * ndim)

    # Grid.dims is the coordinate span along each axis (not the number of cells).
    span = tuple(extent[1][i] - extent[0][i] for i in range(ndim))
    step = tuple(
        span[i] / s if isinstance(s, int) else s
        for i, s in enumerate(size)
    )

    # define grid centerpoint
    kwds['center'] = tuple((extent[0][i] + extent[1][i]) / 2 for i in range(ndim))

    return Grid(span, step, **kwds)

Construct a grid using the provided extent (or the bounds of passed point set) and voxel size (float) or grid dimensions (int).

Parameters

extent : tuple | np.ndarray

The extent of the grid in each dimension as [(xmin,ymin,…), (xmax, ymax, …)]. If a tuple is provided, it is interpreted as the bounds of the grid in each dimension. If a numpy array is provided, the min and max of these points in each dimension will be computed and used as the bounds.

size : tuple

The size of the grid in each dimension. If a float is provided, it is interpreted as the voxel size in each dimension. If an integer is provided, it is interpreted as the number of grid cells in each dimension.

**kwds : dict

Additional keyword arguments are passed directly to Grid.

def poisson_disk_indices_2d(x: numpy.ndarray,
y: numpy.ndarray,
radius: float,
max_points: int,
seed: int | None = None) ‑> numpy.ndarray

Expand source code

def poisson_disk_indices_2d(
    x: np.ndarray,
    y: np.ndarray,
    radius: float,
    max_points: int,
    seed: Optional[int] = None,
    ) -> np.ndarray:
    """
    Fast Poisson-disk sampling on an arbitrary set of (x,y) points.

    A uniform grid-hashing approach ensures that any two chosen points are
    at least ``radius`` apart by checking only the 8 neighboring cells.

    Parameters
    ----------
    x, y : np.ndarray
        1D arrays of the same length N (flattened coordinates).
    radius : float
        Minimum separation between sampled points.
    max_points : int
        Maximum number of points to return (may return fewer if not feasible).
    seed : int | None, default=None
        RNG seed for reproducibility.

    Returns
    -------
    np.ndarray
        Indices (into ``x``/``y``) of selected points, dtype=int64.
    """
    if x.shape != y.shape:
        raise ValueError("x and y must have identical shapes.")
    if radius <= 0:
        raise ValueError("radius must be > 0.")
    if max_points <= 0:
        return np.empty(0, dtype=np.int64)

    rng = np.random.default_rng(seed)
    coords = np.c_[x, y]
    n = coords.shape[0]

    cell = radius / np.sqrt(2.0)  # r/√2 -> 8-neighborhood suffices
    xmin, ymin = coords.min(axis=0)
    ix = np.floor((coords[:, 0] - xmin) / cell).astype(np.int32)
    iy = np.floor((coords[:, 1] - ymin) / cell).astype(np.int32)

    order = rng.permutation(n)
    grid: dict[int, dict[int, int]] = {}
    chosen: List[int] = []
    r2 = radius * radius

    neighbors = [(-1,-1),(-1,0),(-1,1),
                 ( 0,-1),( 0,0),( 0,1),
                 ( 1,-1),( 1,0),( 1,1)]

    for p in order:
        if len(chosen) >= max_points:
            break
        cx, cy = int(ix[p]), int(iy[p])
        ok = True
        for dx, dy in neighbors:
            gx, gy = cx + dx, cy + dy
            col = grid.get(gx)
            if col is None:
                continue
            q = col.get(gy)
            if q is None:
                continue
            dxv = coords[p, 0] - coords[q, 0]
            dyv = coords[p, 1] - coords[q, 1]
            if dxv * dxv + dyv * dyv < r2:
                ok = False
                break
        if ok:
            grid.setdefault(cx, {})[cy] = p
            chosen.append(p)

    return np.asarray(chosen, dtype=np.int64)

Fast Poisson-disk sampling on an arbitrary set of (x,y) points.

A uniform grid-hashing approach ensures that any two chosen points are at least radius apart by checking only the 8 neighboring cells.

Parameters

x, y : np.ndarray

1D arrays of the same length N (flattened coordinates).

radius : float

Minimum separation between sampled points.

max_points : int

Maximum number of points to return (may return fewer if not feasible).

seed : int | None, default=None

RNG seed for reproducibility.

Returns

np.ndarray

Indices (into x/y) of selected points, dtype=int64.

def poisson_disk_indices_3d(x: numpy.ndarray,
y: numpy.ndarray,
z: numpy.ndarray,
radius: float,
max_points: int,
seed: int | None = None) ‑> numpy.ndarray

Expand source code

def poisson_disk_indices_3d(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    radius: float,
    max_points: int,
    seed: Optional[int] = None,
) -> np.ndarray:
    """
    Fast Poisson-disk sampling on an arbitrary set of (x,y,z) points.

    Uses a uniform 3D grid-hashing approach so that any two chosen points are
    at least `radius` apart by checking only the 27 neighboring cells.

    Parameters
    ----------
    x, y, z : np.ndarray
        1D arrays of the same length N (flattened coordinates).
    radius : float
        Minimum separation between sampled points.
    max_points : int
        Maximum number of points to return (may return fewer if not feasible).
    seed : int | None, default=None
        RNG seed for reproducibility.

    Returns
    -------
    np.ndarray
        Indices (into x/y/z) of selected points, dtype=int64.
    """
    if not (x.shape == y.shape == z.shape):
        raise ValueError("x, y, and z must have identical shapes.")
    if radius <= 0:
        raise ValueError("radius must be > 0.")
    if max_points <= 0:
        return np.empty(0, dtype=np.int64)

    rng = np.random.default_rng(seed)
    coords = np.c_[x, y, z]
    n = coords.shape[0]

    # In 3D, choose cube side so the space diagonal equals `radius`:
    # s * sqrt(3) = radius  ->  s = radius / sqrt(3)
    cell = radius / np.sqrt(3.0)

    xmin, ymin, zmin = coords.min(axis=0)
    ix = np.floor((coords[:, 0] - xmin) / cell).astype(np.int32)
    iy = np.floor((coords[:, 1] - ymin) / cell).astype(np.int32)
    iz = np.floor((coords[:, 2] - zmin) / cell).astype(np.int32)

    order = rng.permutation(n)
    grid: dict[Tuple[int, int, int], int] = {}
    chosen: List[int] = []
    r2 = radius * radius

    for p in order:
        if len(chosen) >= max_points:
            break

        cx, cy, cz = int(ix[p]), int(iy[p]), int(iz[p])
        ok = True

        # Check 27 neighboring cells: (dx, dy, dz) in {-1, 0, 1}
        for dx in (-1, 0, 1):
            if not ok:
                break
            for dy in (-1, 0, 1):
                if not ok:
                    break
                for dz in (-1, 0, 1):
                    q = grid.get((cx + dx, cy + dy, cz + dz))
                    if q is None:
                        continue
                    dxv = coords[p, 0] - coords[q, 0]
                    dyv = coords[p, 1] - coords[q, 1]
                    dzv = coords[p, 2] - coords[q, 2]
                    if dxv * dxv + dyv * dyv + dzv * dzv < r2:
                        ok = False
                        break

        if ok:
            grid[(cx, cy, cz)] = p
            chosen.append(p)

    return np.asarray(chosen, dtype=np.int64)

Fast Poisson-disk sampling on an arbitrary set of (x,y,z) points.

Uses a uniform 3D grid-hashing approach so that any two chosen points are at least radius apart by checking only the 27 neighboring cells.

Parameters

x, y, z : np.ndarray

1D arrays of the same length N (flattened coordinates).

radius : float

Minimum separation between sampled points.

max_points : int

Maximum number of points to return (may return fewer if not feasible).

seed : int | None, default=None

RNG seed for reproducibility.

Returns

np.ndarray

Indices (into x/y/z) of selected points, dtype=int64.

def section(dims: tuple,
origin: numpy.ndarray,
normal: numpy.ndarray,
width: float = None,
height: float = None,
step=None)

Expand source code

def section(dims : tuple, origin : np.ndarray, normal : np.ndarray, width : float = None, height : float = None, step=None):
    """
    Create a grid of 3D points on a plane. Useful for evaluating and plotting
    cross sections through 3D models.

    Parameters
    ----------
    dims : tuple of int
        A tuple (nx, ny) specifying the number of grid points in the horizontal and vertical directions.
    origin : array-like, shape (3,)
        The center point of the grid.
    normal : array-like, shape (3,)
        The normal vector defining the plane (need not be unit length).
    width : float
        The total extent of the grid in the horizontal direction. If None, step must be defined.
    height : float
        The total extent of the grid in the vertical direction. If None, step must be defined.
    step : float, optional
        Spacing between grid points. If provided, the grid points will be spaced by this 
        amount in both in-plane directions, and the extents will be determined by `dims` and `step` 
        rather than by the given `width` and `height`.

    Returns
    -------
    grid : np.ndarray
        An array of shape (nx, ny, 3) representing the grid of 3D points.
    """

    # Ensure origin and normal are numpy arrays of type float.
    origin = np.asarray(origin, dtype=float)
    normal = np.asarray(normal, dtype=float)
    
    # Normalize the normal vector.
    n = normal / np.linalg.norm(normal)
    
    # To construct an orthonormal basis for the plane, we need a vector that is
    # not parallel to n. We choose an arbitrary vector; for example, if n is close
    # to [0,0,1] we can choose [1,0,0], otherwise we choose [0,0,1].
    if np.allclose(n, [0, 0, 1]) or np.allclose(n, [0, 0, -1]):
        arbitrary = np.array([0, 1, 0], dtype=float) # height axis is N-S (map)
    else:
        arbitrary = np.array([0, 0, 1], dtype=float) # height axis is up-down (section)
    
    # First in-plane axis: u = n x arbitrary
    u = np.cross(n, arbitrary)
    u /= np.linalg.norm(u)
    
    # Second in-plane axis: v = n x u (this ensures u, v, n is a right-handed frame)
    v = np.cross(n, u)
    v /= np.linalg.norm(v)
    
    nx, ny = dims
    
    if step is None:
        # If step is not provided, we want the grid to exactly span the given width and height.
        # For nx>1, the spacing along the u-direction is width/(nx-1) so that the points
        # go from -width/2 to +width/2.
        u_coords = np.linspace(-width/2, width/2, nx)
        v_coords = np.linspace(-height/2, height/2, ny)
    else:
        # If step is provided, we ignore width and height.
        # We center the grid so that the middle point is at zero.
        if isinstance(step, tuple) or isinstance(step, list):
            sx = step[0]
            sy = step[1]
        else:
            sx = sy = step
        u_coords = -((nx - 1) / 2) * sx + np.arange(nx) * sx
        v_coords = -((ny - 1) / 2) * sy + np.arange(ny) * sy
    
    # Create the grid points.
    grid = np.empty((nx, ny, 3))
    for i, u_val in enumerate(u_coords):
        for j, v_val in enumerate(v_coords):
            grid[i, j] = origin + u_val * u + v_val * v
    
    return grid.shape, grid.reshape((-1,3))

Create a grid of 3D points on a plane. Useful for evaluating and plotting cross sections through 3D models.

Parameters

dims : tuple of int

A tuple (nx, ny) specifying the number of grid points in the horizontal and vertical directions.

origin : array-like, shape (3,)

The center point of the grid.

normal : array-like, shape (3,)

The normal vector defining the plane (need not be unit length).

width : float

The total extent of the grid in the horizontal direction. If None, step must be defined.

height : float

The total extent of the grid in the vertical direction. If None, step must be defined.

step : float, optional

Spacing between grid points. If provided, the grid points will be spaced by this amount in both in-plane directions, and the extents will be determined by dims and step rather than by the given width and height.

Returns

grid : np.ndarray

An array of shape (nx, ny, 3) representing the grid of 3D points.

def triangle_wave(x, A=1, T=6.283185307179586, n_terms=11)

Expand source code

def triangle_wave(x, A=1, T=2*np.pi, n_terms=11):
    """
    Approximates a triangle wave using a Fourier series.

    Parameters:
    x : array-like
        Input values (time or position).
    A : float
        Amplitude of the triangle wave.
    T : float
        Period (wavelength) of the triangle wave.
    n_terms : int
        Number of terms in the Fourier series approximation.

    Returns:
    f_x : array-like
        Approximated triangle wave values.
    """
    f_x = x * 0
    isTorch = isinstance(f_x, torch.Tensor)
    for n in range(1, n_terms + 1, 2):  # Only odd harmonics
        if isTorch:
            f_x = f_x + ((-1) ** ((n-1)//2)*(A/n**2))*torch.sin((2*torch.pi*n/T)*x)
        else:
            f_x = f_x + ((-1) ** ((n-1)//2)*(A/n**2))*np.sin((2*np.pi*n/T)*x)
    return (8/np.pi**2)*f_x  # Normalizing factor for amplitude scaling

Approximates a triangle wave using a Fourier series.

Parameters: x : array-like Input values (time or position). A : float Amplitude of the triangle wave. T : float Period (wavelength) of the triangle wave. n_terms : int Number of terms in the Fourier series approximation.

Returns: f_x : array-like Approximated triangle wave values.

Classes

class Grid (dims, step, center=None, rotation=None, sampleArgs={})

Expand source code

class Grid(object):
    """
    Class encapsulating more advanced grids (e.g., non-axis aligned grids) and facilitating array reshaping.
    """
    def __init__( self, dims, step, center=None, rotation=None, sampleArgs={} ):
        """
        Initialize a grid object with specified dimensions, step size, origin, and rotation.
        Parameters:
            dims (tuple): Dimensions of the grid. A tuple of integers specifying the size along each axis.
            step (float or int or tuple): Step size for the grid. If a single float or int is provided, 
                it is applied uniformly to all dimensions. If a tuple is provided, it specifies the step 
                size for each dimension.
            center (tuple, optional): The origin (center) of the grid in the global coordinate space. Defaults to the origin.
            rotation (numpy.ndarray, optional): A rotation matrix defining the orientation of the grid. 
                Defaults to None.
            sampleArgs (dict, optional): Additional arguments for sampling with draw(). Defaults to an empty dictionary.
        Attributes:
            dims (tuple): Stored dimensions of the grid.
            step (tuple): Step size for each dimension.
            ndim (int): Number of dimensions in the grid.
            sampleArgs (dict): Additional arguments for sampling.
            axes (list): List of numpy arrays representing the axes of the grid.
            shape (tuple): Shape of the grid based on the axes.
            matrix (numpy.ndarray): Transformation matrix combining origin and rotation.
            center (numpy.ndarray): Stored origin (center) of the grid for convenience.
        """

        self.dims = tuple([i for i in dims]) # store dims
        if isinstance(step, float) or isinstance(step, int):  # expand to step each dimension if a float or int is passed
            step = tuple([step for i in dims])
        self.step = step
        self.ndim = len(self.dims) # number of dimensions in grid (2- or 3)
        self.sampleArgs = sampleArgs

        # build axes arrays
        self.axes = [np.arange(d, step=s) - d/2 for d,s in zip(self.dims, self.step)]
        self.shape = tuple([len(a) for a in self.axes]) # store grid shape
        
        # build matrix
        self.matrix = np.eye(self.ndim+1)
        if center is not None: # define grid origin
            self.matrix[:self.ndim, self.ndim] = np.array(center)
            self.tformed_axes = [self.axes[i] + center[i] for i in range(len(center))]
        if rotation is not None: # define grid rotation matrix
            self.matrix[:self.ndim, :self.ndim] = rotation
        self.center = self.matrix[:self.ndim, self.ndim] # store for convenience
        self._clearCache()

    def __str__(self):
        from curlew.text import grid_str
        return grid_str(self)

    def coords(self, transform=True):
        """
        Get an (N,2) or (N,3) array of coordinates representing the position of each cell in this grid. 
        """
        if self._cache is not None:
            assert transform, "Transform cannot be False if altCoords have been defined."
            return self._cache # easy!

        # create ravelled list of points
        #points = np.array( np.meshgrid( *self.axes, indexing='xy' ) ).T.reshape((-1,self.ndim))
        #points = np.array( np.meshgrid( *self.axes, indexing='xy' ) ).reshape((-1,self.ndim))
        coords = np.meshgrid(*self.axes, indexing='ij')
        grid = np.array([c.T for c in coords[::-1]]).T
        grid = grid[...,::-1]
        points = grid.reshape((-1, self.ndim))

        # apply transform matrix
        if transform:
            points = np.hstack([points, np.ones( (len(points),1) )]) # concatentate additional dimension
            points = (points @ self.matrix.T)[:,:-1] # drop fake dimension
        return points
    
    def reshape(self, values):
        """
        Reshape an (N,) or (N,d) array of points to an N-d grid (image or voxel) array.
        """
        if len(values.shape) == 1:
            return values.reshape(self.shape) # easy
        else:
            return values.reshape( self.shape + (values.shape[-1], ) )
    
    def _setCache(self, coords):
        """
        Store a cached transformed grid. Used during training of individual fields, to avoid lots of recomputations 
        of transformed grid coords.
        """
        assert (np.array( coords.shape ) == np.array( self.coords().shape )).all()
        self._cache = coords
    
    def _clearCache(self):
        """ Remove any defined cache. """
        self._cache = None
    
    def draw(self, transform=None):
        """
        Draw a random sample according to `self.sampleArgs`. This is used to draw
        random samples (using `self.sample`) when fitting global constraints.
        """
        coords = self.sample(**self.sampleArgs, tensor=True)
        if transform and (self._cache is not None): # we need to apply a transform?
            # N.B. assumes that if self._cache is defined then this contains
            # pre-computed transformed grid point positions
            coords = transform(coords, end=None)
        return coords

    def sample(self, N=4096, poissonDisk=None, tensor=False):
        """
        Sample random points from this grid. Used for e.g., evaluating "global" constraints. 

        Parameters
        ----------
        N : int
            The number of random samples to select.
        poissonDisk : tuple | optional (defaults to None)
            If not None, should be given as a tuple (r, k, seed), where r is the minimum distance between
            the points and k is the maximum number of points allowed, and seed ensures reproducibility
            (can be None to use previous random seed).
            If None, draws randomly from the grid.
        tensor : bool | optional
            True if sampled points should be returned as a torch.Tensor. Default is False. 
        """

        # get grid of points to sample from
        grid = self.coords()

        # do random sampling
        if poissonDisk is not None: # Do poisson sampling to ensure evenly spaced points
            if grid.shape[1] == 2: # 2D Grid
                ix = poisson_disk_indices_2d(grid[:, 0], grid[:, 1],
                                              radius=poissonDisk[0], max_points=poissonDisk[1],
                                              seed=poissonDisk[2])
            else: # 3D Grid
                ix = poisson_disk_indices_3d(grid[:, 0], grid[:, 1], grid[:, 2],
                                              radius=poissonDisk[0], max_points=poissonDisk[1],
                                              seed=poissonDisk[2])

        else:
            ix = np.random.choice(len(grid), N, replace=False ) # draw random points from the grid (without replacement)
        
        # get chosen points
        out = grid[ix,:]

        # return in desired format
        if tensor:
            from curlew.core import _tensor
            return _tensor(out)
        else:
            return out
    
    def contour(self, values, iso, normals=False, transform=True, mask=None, erodeMask=0 ):
      """
      Use values computed for this grid to extract 2D (lines) or 3D (surfaces) contours. Requires scikit-image.

      Parameters
      -----------
          values : np.ndarray
              A 1D array of values to extract contours from. Length must be the same as `self.coords()`.
          iso : float
              The isovalue to extract. Must be within the range of `values`.
          normals : bool
              If True, an array of vertex normals will be returned in addition to the vertices and faces.
          transform : bool | optional
              True (default) if this grids transorm matrix should be applied to return coordinates in world coords. Otherwise
              internal (local) grid coordinates will be returned.
          mask : np.ndarray | None
              Optional boolean mask specifying which grid cells are considered valid for contouring.
              Can be provided either as a ravelled 1D array with the same length as `values`, or
              as an array already shaped like this grid (`self.shape`). Masked-out regions are
              ignored by the contouring algorithm.
          erodeMask : int
              If non-zero and `mask` is provided, morphologically filter the combined mask
              before contouring. Positive values apply binary erosion (trim boundary voxels);
              negative values apply binary dilation (expand the mask by ``abs(erodeMask)``
              iterations). Dilation is often needed when the mask boundary coincides with the
              isosurface being extracted. Default is 0 (no filtering).
      
      Returns
      --------
          An array of `vert`, `faces` for 3D grids, or a list of contour polylines (lists of coordinates) for 2D grids.
      """
      try:
          import skimage
      except:
          assert False, "Please install scikit-image using `pip install scikit-image`"
      
      vol = self.reshape(values)
      contour_mask = np.isfinite(vol)
      if mask is not None:
          m = np.asarray(mask)
          if m.dtype != bool:
              m = m.astype(bool)
          if m.shape == (len(values),):
              mask = self.reshape(m)
          elif m.shape == self.shape:
              mask = m
          else:
              raise ValueError(
                  "volume_mask must be either 1D (len(values),) or shaped like the grid (self.shape). "
                  f"Got volume_mask.shape={m.shape}, expected {(len(values),)} or {self.shape}."
              )
          contour_mask = contour_mask & mask
          if erodeMask:
              if erodeMask > 0:
                  from scipy.ndimage import binary_erosion
                  contour_mask = binary_erosion(contour_mask, iterations=int(erodeMask))
              else:
                  from scipy.ndimage import binary_dilation
                  contour_mask = binary_dilation(contour_mask, iterations=int(-erodeMask))

      if self.ndim == 3: # marching cubes
          from skimage.measure import marching_cubes
          vix, faces, norm, _ = marching_cubes( vol, level=iso, mask=contour_mask ) # find isosurface using marching cubes 

          # build interpolators that convert indices to positions
          from scipy.interpolate import interp1d
          ix = [ interp1d(np.arange(self.shape[i]), self.axes[i]) for i in range(len(self.shape)) ]
          verts = np.vstack( [ix[i]( vix[:, i] ) for i in range(len(self.shape))] ).T

          # apply rotation matrix
          if transform:
              verts = np.hstack([verts, np.ones( (len(verts),1) )]) # concatentate fourth dimension
              verts = (verts @ self.matrix.T)[:,:-1] # drop fake dimension

          # remove invalid vertices
          if np.isnan(verts).any():
              vmask = ~np.isnan(verts).any(axis=-1)
              verts = verts[vmask]
              norm = norm[vmask]
                
              # encode faces again...
              old2New=np.arange(len(vmask), dtype=float)
              old2New[vmask] = np.arange(np.sum(vmask))
              old2New[~vmask] = np.nan
              faces = np.array( [(old2New[i], old2New[j], old2New[k]) for i,j,k in faces] )
              faces = faces[np.isfinite(faces).all(axis=-1)]

          if normals:
              return verts, faces, norm
          else:
              return verts, faces
      elif self.ndim == 2: # marching squares
          from skimage.measure import find_contours
          contours = find_contours( vol, iso, mask=contour_mask)

          # build interpolators that convert indices to positions
          from scipy.interpolate import interp1d
          for i,vix in enumerate(contours):
              ix = [ interp1d(np.arange(self.shape[i]), self.axes[i], fill_value='extrapolate') for i in range(len(self.shape)) ]
              contours[i] = np.vstack( [ix[i]( vix[:, i] ) for i in range(len(self.shape))] ).T

          if transform:
              for i,c in enumerate(contours):
                points = np.hstack([c, np.ones( (len(c),1) )]) # concatentate additional dimension
                contours[i] = (points @ self.matrix.T)[:,:-1] # apply and drop fake dimension

          # TODO; mask NaN areas from 2D contours

          return contours
        
    def copy(self):
        """
        Create a copy of the current Grid instance.
        """
        new_grid = Grid(self.dims, self.step, center=self.center, rotation=self.matrix[:self.ndim, :self.ndim], sampleArgs=self.sampleArgs)
        if self._cache is not None:
            new_grid._setCache(self._cache)
        return new_grid

Class encapsulating more advanced grids (e.g., non-axis aligned grids) and facilitating array reshaping.

Initialize a grid object with specified dimensions, step size, origin, and rotation.

Parameters

dims (tuple): Dimensions of the grid. A tuple of integers specifying the size along each axis. step (float or int or tuple): Step size for the grid. If a single float or int is provided, it is applied uniformly to all dimensions. If a tuple is provided, it specifies the step size for each dimension. center (tuple, optional): The origin (center) of the grid in the global coordinate space. Defaults to the origin. rotation (numpy.ndarray, optional): A rotation matrix defining the orientation of the grid. Defaults to None. sampleArgs (dict, optional): Additional arguments for sampling with draw(). Defaults to an empty dictionary.

Attributes

dims : tuple

Stored dimensions of the grid.

step : tuple

Step size for each dimension.

ndim : int

Number of dimensions in the grid.

sampleArgs : dict

Additional arguments for sampling.

axes : list

List of numpy arrays representing the axes of the grid.

shape : tuple

Shape of the grid based on the axes.

matrix : numpy.ndarray

Transformation matrix combining origin and rotation.

center : numpy.ndarray

Stored origin (center) of the grid for convenience.

Methods

def contour(self, values, iso, normals=False, transform=True, mask=None, erodeMask=0)

Expand source code

def contour(self, values, iso, normals=False, transform=True, mask=None, erodeMask=0 ):
  """
  Use values computed for this grid to extract 2D (lines) or 3D (surfaces) contours. Requires scikit-image.

  Parameters
  -----------
      values : np.ndarray
          A 1D array of values to extract contours from. Length must be the same as `self.coords()`.
      iso : float
          The isovalue to extract. Must be within the range of `values`.
      normals : bool
          If True, an array of vertex normals will be returned in addition to the vertices and faces.
      transform : bool | optional
          True (default) if this grids transorm matrix should be applied to return coordinates in world coords. Otherwise
          internal (local) grid coordinates will be returned.
      mask : np.ndarray | None
          Optional boolean mask specifying which grid cells are considered valid for contouring.
          Can be provided either as a ravelled 1D array with the same length as `values`, or
          as an array already shaped like this grid (`self.shape`). Masked-out regions are
          ignored by the contouring algorithm.
      erodeMask : int
          If non-zero and `mask` is provided, morphologically filter the combined mask
          before contouring. Positive values apply binary erosion (trim boundary voxels);
          negative values apply binary dilation (expand the mask by ``abs(erodeMask)``
          iterations). Dilation is often needed when the mask boundary coincides with the
          isosurface being extracted. Default is 0 (no filtering).
  
  Returns
  --------
      An array of `vert`, `faces` for 3D grids, or a list of contour polylines (lists of coordinates) for 2D grids.
  """
  try:
      import skimage
  except:
      assert False, "Please install scikit-image using `pip install scikit-image`"
  
  vol = self.reshape(values)
  contour_mask = np.isfinite(vol)
  if mask is not None:
      m = np.asarray(mask)
      if m.dtype != bool:
          m = m.astype(bool)
      if m.shape == (len(values),):
          mask = self.reshape(m)
      elif m.shape == self.shape:
          mask = m
      else:
          raise ValueError(
              "volume_mask must be either 1D (len(values),) or shaped like the grid (self.shape). "
              f"Got volume_mask.shape={m.shape}, expected {(len(values),)} or {self.shape}."
          )
      contour_mask = contour_mask & mask
      if erodeMask:
          if erodeMask > 0:
              from scipy.ndimage import binary_erosion
              contour_mask = binary_erosion(contour_mask, iterations=int(erodeMask))
          else:
              from scipy.ndimage import binary_dilation
              contour_mask = binary_dilation(contour_mask, iterations=int(-erodeMask))

  if self.ndim == 3: # marching cubes
      from skimage.measure import marching_cubes
      vix, faces, norm, _ = marching_cubes( vol, level=iso, mask=contour_mask ) # find isosurface using marching cubes 

      # build interpolators that convert indices to positions
      from scipy.interpolate import interp1d
      ix = [ interp1d(np.arange(self.shape[i]), self.axes[i]) for i in range(len(self.shape)) ]
      verts = np.vstack( [ix[i]( vix[:, i] ) for i in range(len(self.shape))] ).T

      # apply rotation matrix
      if transform:
          verts = np.hstack([verts, np.ones( (len(verts),1) )]) # concatentate fourth dimension
          verts = (verts @ self.matrix.T)[:,:-1] # drop fake dimension

      # remove invalid vertices
      if np.isnan(verts).any():
          vmask = ~np.isnan(verts).any(axis=-1)
          verts = verts[vmask]
          norm = norm[vmask]
            
          # encode faces again...
          old2New=np.arange(len(vmask), dtype=float)
          old2New[vmask] = np.arange(np.sum(vmask))
          old2New[~vmask] = np.nan
          faces = np.array( [(old2New[i], old2New[j], old2New[k]) for i,j,k in faces] )
          faces = faces[np.isfinite(faces).all(axis=-1)]

      if normals:
          return verts, faces, norm
      else:
          return verts, faces
  elif self.ndim == 2: # marching squares
      from skimage.measure import find_contours
      contours = find_contours( vol, iso, mask=contour_mask)

      # build interpolators that convert indices to positions
      from scipy.interpolate import interp1d
      for i,vix in enumerate(contours):
          ix = [ interp1d(np.arange(self.shape[i]), self.axes[i], fill_value='extrapolate') for i in range(len(self.shape)) ]
          contours[i] = np.vstack( [ix[i]( vix[:, i] ) for i in range(len(self.shape))] ).T

      if transform:
          for i,c in enumerate(contours):
            points = np.hstack([c, np.ones( (len(c),1) )]) # concatentate additional dimension
            contours[i] = (points @ self.matrix.T)[:,:-1] # apply and drop fake dimension

      # TODO; mask NaN areas from 2D contours

      return contours

Use values computed for this grid to extract 2D (lines) or 3D (surfaces) contours. Requires scikit-image.

Parameters

values : np.ndarray
    A 1D array of values to extract contours from. Length must be the same as <code>self.coords()</code>.
iso : float
    The isovalue to extract. Must be within the range of <code>values</code>.
normals : bool
    If True, an array of vertex normals will be returned in addition to the vertices and faces.
transform : bool | optional
    True (default) if this grids transorm matrix should be applied to return coordinates in world coords. Otherwise
    internal (local) grid coordinates will be returned.
mask : np.ndarray | None
    Optional boolean mask specifying which grid cells are considered valid for contouring.
    Can be provided either as a ravelled 1D array with the same length as <code>values</code>, or
    as an array already shaped like this grid (<code>self.shape</code>). Masked-out regions are
    ignored by the contouring algorithm.
erodeMask : int
    If non-zero and <code>mask</code> is provided, morphologically filter the combined mask
    before contouring. Positive values apply binary erosion (trim boundary voxels);
    negative values apply binary dilation (expand the mask by <code>abs(erodeMask)</code>
    iterations). Dilation is often needed when the mask boundary coincides with the
    isosurface being extracted. Default is 0 (no filtering).

Returns

An array of <code>vert</code>, <code>faces</code> for 3D grids, or a list of contour polylines (lists of coordinates) for 2D grids.

def coords(self, transform=True)

Expand source code

def coords(self, transform=True):
    """
    Get an (N,2) or (N,3) array of coordinates representing the position of each cell in this grid. 
    """
    if self._cache is not None:
        assert transform, "Transform cannot be False if altCoords have been defined."
        return self._cache # easy!

    # create ravelled list of points
    #points = np.array( np.meshgrid( *self.axes, indexing='xy' ) ).T.reshape((-1,self.ndim))
    #points = np.array( np.meshgrid( *self.axes, indexing='xy' ) ).reshape((-1,self.ndim))
    coords = np.meshgrid(*self.axes, indexing='ij')
    grid = np.array([c.T for c in coords[::-1]]).T
    grid = grid[...,::-1]
    points = grid.reshape((-1, self.ndim))

    # apply transform matrix
    if transform:
        points = np.hstack([points, np.ones( (len(points),1) )]) # concatentate additional dimension
        points = (points @ self.matrix.T)[:,:-1] # drop fake dimension
    return points

Get an (N,2) or (N,3) array of coordinates representing the position of each cell in this grid.

def copy(self)

Expand source code

def copy(self):
    """
    Create a copy of the current Grid instance.
    """
    new_grid = Grid(self.dims, self.step, center=self.center, rotation=self.matrix[:self.ndim, :self.ndim], sampleArgs=self.sampleArgs)
    if self._cache is not None:
        new_grid._setCache(self._cache)
    return new_grid

Create a copy of the current Grid instance.

def draw(self, transform=None)

Expand source code

def draw(self, transform=None):
    """
    Draw a random sample according to `self.sampleArgs`. This is used to draw
    random samples (using `self.sample`) when fitting global constraints.
    """
    coords = self.sample(**self.sampleArgs, tensor=True)
    if transform and (self._cache is not None): # we need to apply a transform?
        # N.B. assumes that if self._cache is defined then this contains
        # pre-computed transformed grid point positions
        coords = transform(coords, end=None)
    return coords

Draw a random sample according to self.sampleArgs. This is used to draw random samples (using self.sample) when fitting global constraints.

def reshape(self, values)

Expand source code

def reshape(self, values):
    """
    Reshape an (N,) or (N,d) array of points to an N-d grid (image or voxel) array.
    """
    if len(values.shape) == 1:
        return values.reshape(self.shape) # easy
    else:
        return values.reshape( self.shape + (values.shape[-1], ) )

Reshape an (N,) or (N,d) array of points to an N-d grid (image or voxel) array.

def sample(self, N=4096, poissonDisk=None, tensor=False)

Expand source code

def sample(self, N=4096, poissonDisk=None, tensor=False):
    """
    Sample random points from this grid. Used for e.g., evaluating "global" constraints. 

    Parameters
    ----------
    N : int
        The number of random samples to select.
    poissonDisk : tuple | optional (defaults to None)
        If not None, should be given as a tuple (r, k, seed), where r is the minimum distance between
        the points and k is the maximum number of points allowed, and seed ensures reproducibility
        (can be None to use previous random seed).
        If None, draws randomly from the grid.
    tensor : bool | optional
        True if sampled points should be returned as a torch.Tensor. Default is False. 
    """

    # get grid of points to sample from
    grid = self.coords()

    # do random sampling
    if poissonDisk is not None: # Do poisson sampling to ensure evenly spaced points
        if grid.shape[1] == 2: # 2D Grid
            ix = poisson_disk_indices_2d(grid[:, 0], grid[:, 1],
                                          radius=poissonDisk[0], max_points=poissonDisk[1],
                                          seed=poissonDisk[2])
        else: # 3D Grid
            ix = poisson_disk_indices_3d(grid[:, 0], grid[:, 1], grid[:, 2],
                                          radius=poissonDisk[0], max_points=poissonDisk[1],
                                          seed=poissonDisk[2])

    else:
        ix = np.random.choice(len(grid), N, replace=False ) # draw random points from the grid (without replacement)
    
    # get chosen points
    out = grid[ix,:]

    # return in desired format
    if tensor:
        from curlew.core import _tensor
        return _tensor(out)
    else:
        return out

Sample random points from this grid. Used for e.g., evaluating "global" constraints.

Parameters

N : int

The number of random samples to select.

poissonDisk : tuple | optional (defaults to None)

If not None, should be given as a tuple (r, k, seed), where r is the minimum distance between the points and k is the maximum number of points allowed, and seed ensures reproducibility (can be None to use previous random seed). If None, draws randomly from the grid.

tensor : bool | optional

True if sampled points should be returned as a torch.Tensor. Default is False.

class Transform (matrix: numpy.ndarray | ForwardRef('torch.Tensor'))

Expand source code

@dataclass
class Transform:
    """
    Homogeneous transform supporting:
      - 3x3 matrices for 2D points
      - 4x4 matrices for 3D points
    """
    matrix: ArrayLike

    def __post_init__(self):
        if isinstance(self.matrix, int):
            self.matrix = np.eye(self.matrix + 1)

        # If handed a CUDA tensor, pull it to CPU before storing as numpy
        if isinstance(self.matrix, torch.Tensor):
            self.matrix = self.matrix.detach().cpu().numpy()

        self.matrix = _numpy(self.matrix)

        if self.matrix.shape not in [(3, 3), (4, 4)]:
            raise ValueError("Transform matrix must be 3x3 or 4x4")

    def isIdentity(self) -> bool:
        """Return True if this transform is the identity matrix."""
        return np.allclose(self.matrix, np.eye(self.matrix.shape[0]))

    def set(self, mat: ArrayLike):
        """Set the transform matrix."""
        self.matrix = mat
        self.__post_init__()

    def apply(self, points: ArrayLike) -> ArrayLike:
        """
        Apply the transform to an (N,2) or (N,3) array of points.
        Handles both numpy arrays and pytorch tensors.
        """
        if isinstance(points, torch.Tensor):
            return self._apply_torch(points)
        elif isinstance(points, np.ndarray):
            return self._apply_numpy(points)
        else:
            raise TypeError("points must be a numpy array or torch tensor")

    def __call__(self, points: ArrayLike) -> ArrayLike:
        """Alias for apply()."""
        return self.apply(points)

    def inverse(self) -> "Transform":
        """Return a Transform representing the inverse operation."""
        # Always invert on CPU to avoid cuSOLVER initialisation issues;
        # _apply_torch will move the matrix to the right device on demand.
        mat = torch.as_tensor(self.matrix)  # CPU, from numpy source of truth

        try:
            inv = torch.linalg.inv(mat)
        except torch.linalg.LinAlgError as e:
            raise ValueError("Transform matrix is not invertible") from e
        # Any other RuntimeError (CUDA, dtype, driver) propagates naturally

        return Transform(inv.numpy())

    def _apply_numpy(self, points: np.ndarray) -> np.ndarray:
        """Apply the transform to a numpy (N, D) array."""
        if points.ndim != 2:
            raise ValueError("points must have shape (N, D)")

        n, d = points.shape
        mat = self.matrix

        if mat.shape == (3, 3) and d != 2:
            raise ValueError("3x3 transform requires (N,2) points")
        if mat.shape == (4, 4) and d != 3:
            raise ValueError("4x4 transform requires (N,3) points")

        ones = np.ones((n, 1), dtype=points.dtype)
        homo = np.concatenate([points, ones], axis=1)   # (N, D+1)
        transformed = (mat @ homo.T).T                   # (N, D+1)
        w = transformed[:, -1:]
        return transformed[:, :-1] / w

    def _apply_torch(self, points: torch.Tensor) -> torch.Tensor:
        """Apply the transform to a torch (N, D) tensor."""
        if points.ndim != 2:
            raise ValueError("points must have shape (N, D)")

        n, d = points.shape

        # Cast lazily: matches device and dtype of the incoming points.
        # Zero-copy if already matching; cheap cast otherwise.
        # Keeps self.matrix as the CPU numpy source of truth.
        mat = torch.as_tensor(self.matrix, dtype=points.dtype, device=points.device)

        if mat.shape == (3, 3) and d != 2:
            raise ValueError("3x3 transform requires (N,2) points")
        if mat.shape == (4, 4) and d != 3:
            raise ValueError("4x4 transform requires (N,3) points")

        # Mirror device and dtype of points exactly
        ones = torch.ones((n, 1), dtype=points.dtype, device=points.device)
        homo = torch.cat([points, ones], dim=1)           # (N, D+1)
        transformed = (mat @ homo.T).T                    # (N, D+1)
        w = transformed[:, -1:].clamp(min=1e-12)
        return transformed[:, :-1] / w

Homogeneous transform supporting: - 3x3 matrices for 2D points - 4x4 matrices for 3D points

Instance variables

var matrix : numpy.ndarray | torch.Tensor

Methods

def apply(self, points: numpy.ndarray | ForwardRef('torch.Tensor')) ‑> numpy.ndarray | torch.Tensor

Expand source code

def apply(self, points: ArrayLike) -> ArrayLike:
    """
    Apply the transform to an (N,2) or (N,3) array of points.
    Handles both numpy arrays and pytorch tensors.
    """
    if isinstance(points, torch.Tensor):
        return self._apply_torch(points)
    elif isinstance(points, np.ndarray):
        return self._apply_numpy(points)
    else:
        raise TypeError("points must be a numpy array or torch tensor")

Apply the transform to an (N,2) or (N,3) array of points. Handles both numpy arrays and pytorch tensors.

def inverse(self) ‑> Transform

Expand source code

def inverse(self) -> "Transform":
    """Return a Transform representing the inverse operation."""
    # Always invert on CPU to avoid cuSOLVER initialisation issues;
    # _apply_torch will move the matrix to the right device on demand.
    mat = torch.as_tensor(self.matrix)  # CPU, from numpy source of truth

    try:
        inv = torch.linalg.inv(mat)
    except torch.linalg.LinAlgError as e:
        raise ValueError("Transform matrix is not invertible") from e
    # Any other RuntimeError (CUDA, dtype, driver) propagates naturally

    return Transform(inv.numpy())

Return a Transform representing the inverse operation.

def isIdentity(self) ‑> bool

Expand source code

def isIdentity(self) -> bool:
    """Return True if this transform is the identity matrix."""
    return np.allclose(self.matrix, np.eye(self.matrix.shape[0]))

Return True if this transform is the identity matrix.

def set(self, mat: numpy.ndarray | ForwardRef('torch.Tensor'))

Expand source code

def set(self, mat: ArrayLike):
    """Set the transform matrix."""
    self.matrix = mat
    self.__post_init__()

Set the transform matrix.