Skip to content

grid

Module documentation for grid.

API Reference

average_grids(grids, mask='nan')

Compute the average of multiple grids.

:param grids: The input grids to average. :type grids: list of ndarrays :param mask: The mask type to use for averaging. If 'nan', only non-NaN values are considered. Default is 'nan'. :type mask: str :return: The average grid. :rtype: ndarray :raises NotImplementedError: If the mask option is not implemented.

Example

grid1 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) grid2 = np.array([[2, 4, 6], [8, 10, 12], [14, 16, 18]]) averaged_grid = average_grids([grid1, grid2], mask='nan') print(averaged_grid)

Source code in mrdna/arbdmodel/submodule/grid.py
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def average_grids(grids, mask='nan'):
    """
    Compute the average of multiple grids.

    :param grids: The input grids to average.
    :type grids: list of ndarrays
    :param mask: The mask type to use for averaging. If 'nan', only non-NaN values are considered. Default is 'nan'.
    :type mask: str
    :return: The average grid.
    :rtype: ndarray
    :raises NotImplementedError: If the mask option is not implemented.

    Example:
        >>> grid1 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
        >>> grid2 = np.array([[2, 4, 6], [8, 10, 12], [14, 16, 18]])
        >>> averaged_grid = average_grids([grid1, grid2], mask='nan')
        >>> print(averaged_grid)
    """

    def __get_mask(grid, mask):
      if mask == 'nan':
        g_mask = ~np.isnan(grid)
      else:
        raise NotImplementedError
      return g_mask

    counts = np.zeros(grids[0].shape, dtype=int)
    totals = np.zeros(counts.shape)
    for g in grids:
      g_mask = __get_mask(g, mask)
      counts = counts + g_mask
      totals[g_mask] = totals[g_mask] + g[g_mask]
    sl = counts > 0

    average = totals
    average[sl] = totals[sl] / counts[sl]
    average[counts == 0] = np.nan
    return average

convolve_kernel_truncate(array, kernel)

Convolve an array with a kernel, truncating the output.

Parameters:

Name Type Description Default
array ndarray

The input array.

required
kernel ndarray

The kernel to convolve with the array.

required

Returns:

Name Type Description
ndarray

The truncated convolution result.

Raises:

Type Description
AssertionError

If the dimensions of the array and kernel do not match.

AssertionError

If any dimension of the array is smaller than the corresponding dimension of the kernel.

Notes
  • This function assumes that the kernel has odd dimensions along each dimension.
  • If any dimension of the kernel has an even number of elements, a warning message is printed to stderr indicating that the output may be shifted.
Example

array = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) kernel = np.array([[0.5, 0.5], [0.5, 0.5]]) result = convolve_kernel_truncate(array, kernel) print(result)

Source code in mrdna/arbdmodel/submodule/grid.py
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def convolve_kernel_truncate( array, kernel ):
    """
    Convolve an array with a kernel, truncating the output.

    Parameters:
        array (ndarray): The input array.
        kernel (ndarray): The kernel to convolve with the array.

    Returns:
        ndarray: The truncated convolution result.

    Raises:
        AssertionError: If the dimensions of the array and kernel do not match.
        AssertionError: If any dimension of the array is smaller than the corresponding dimension of the kernel.

    Notes:
        - This function assumes that the kernel has odd dimensions along each dimension.
        - If any dimension of the kernel has an even number of elements, a warning message is printed to stderr indicating that the output may be shifted.

    Example:
        >>> array = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
        >>> kernel = np.array([[0.5, 0.5], [0.5, 0.5]])
        >>> result = convolve_kernel_truncate(array, kernel)
        >>> print(result)
    """
    arrayShape = np.shape( array )
    kernelShape = np.shape( kernel )
    assert( len(arrayShape) == len(kernelShape) )

    for an,kn in zip(arrayShape,kernelShape): assert(an > kn)

    dim = 0
    for an in kernelShape:
        dim += 1
        if an % 2 == 0: 
            print("WARNING: kernel has even number of elements along dimension %d; Output may be shifted." % dim, file=sys.stderr)



    initShape = [a+b for a,b in zip(arrayShape,kernelShape)]
    convolution = np.zeros(initShape)
    count = np.zeros(initShape)

    for idx in range( np.prod(kernelShape) ):

        ijk = [ int(idx/a) % b for a,b in
                  zip( np.hstack(([1],np.cumprod(kernelShape[:-1]))), kernelShape) ][::-1]

        s = tuple([ slice( (kn-i), an+kn-i ) for i,kn,an in 
                    zip(ijk,kernelShape,arrayShape) ])

        val = kernel.flatten()[idx]

        convolution[s] += val*array
        count[s] += 1

    ids = count > 0
    convolution[ids] = (convolution[ids]/count[ids]) * np.prod(kernelShape)


    s = tuple([ slice( int((kn-1)/2)+1, -int((kn-1)/2) ) for kn in kernelShape ])
    return convolution[s]

create_bounding_grid(*grids)

Construct a grid bigger than all the (GridDataFormats) grids provided in the arguments; note that the inputs must be orthonormal and have exactly overlapping voxels

Source code in mrdna/arbdmodel/submodule/grid.py
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
def create_bounding_grid( *grids ):
    """ Construct a grid bigger than all the (GridDataFormats) grids provided in the arguments; note that the inputs must be orthonormal and have exactly overlapping voxels """
    def _create_bounding_grid( grid1, grid2 ):
      assert( len(grid1.grid.shape) == 3 )
      assert( len(grid2.grid.shape) == 3 )
      assert(len(grid1.delta) == 3)
      assert( np.all(grid1.delta == grid2.delta) )
      min_origin = np.minimum(*[g.origin for g in (grid1, grid2)])
      max_shape = np.maximum(*[(g.origin-min_origin)//g.delta + np.array(g.grid.shape) for g in (grid1, grid2)]).astype(int)
      assert( np.all( max_shape > 0) )

      new_grid_data = np.empty(max_shape)
      new_grid = Grid( new_grid_data, origin = min_origin, delta = grid1.delta )
      return new_grid

    g1 = grids[0]
    if len( grids ) == 1:
      g1 = _create_bounding_grid( g1, g1 )
    else:
      for g2 in grids[1:]:
        g1 = _create_bounding_grid( g1, g2 )
    g1.grid[...,:] = 0

    return g1

fill_nans(grid, neighborhood=1, max_iterations=np.inf, mask=None)

Fill NaN values in a grid using neighborhood averaging.

Parameters:

Name Type Description Default
grid ndarray

The input grid containing NaN values.

required
neighborhood int

The size of the neighborhood to use for averaging. Default is 1.

1
max_iterations int

The maximum number of iterations to perform. Default is infinity.

inf
mask ndarray

A mask specifying which values to consider for filling NaNs. If None, all non-NaN values are used. Default is None.

None

Returns:

Name Type Description
ndarray

The grid with NaN values filled using neighborhood averaging.

Notes
  • This function requires the 'skimage' package, specifically the 'find_boundaries' function from 'skimage.segmentation'.
Example

grid = np.array([[1, 2, np.nan], [4, np.nan, 6], [np.nan, 8, 9]]) filled_grid = fill_nans(grid, neighborhood=2, max_iterations=10) print(filled_grid)

Source code in mrdna/arbdmodel/submodule/grid.py
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
def fill_nans(grid, neighborhood=1, max_iterations=np.inf, mask=None):
    """
    Fill NaN values in a grid using neighborhood averaging.

    Parameters:
        grid (ndarray): The input grid containing NaN values.
        neighborhood (int, optional): The size of the neighborhood to use for averaging. Default is 1.
        max_iterations (int, optional): The maximum number of iterations to perform. Default is infinity.
        mask (ndarray, optional): A mask specifying which values to consider for filling NaNs. If None, all non-NaN values are used. Default is None.

    Returns:
        ndarray: The grid with NaN values filled using neighborhood averaging.

    Notes:
        - This function requires the 'skimage' package, specifically the 'find_boundaries' function from 'skimage.segmentation'.

    Example:
        >>> grid = np.array([[1, 2, np.nan], [4, np.nan, 6], [np.nan, 8, 9]])
        >>> filled_grid = fill_nans(grid, neighborhood=2, max_iterations=10)
        >>> print(filled_grid)
    """
    from skimage.segmentation import find_boundaries

    ret = np.array(grid)
    if mask is None: mask = ~np.isnan(grid)
    assert(np.all(np.isnan(mask)) == 0)
    nans= np.isnan(ret)

    i = 0
    while np.sum(nans) > 0:
      if i > max_iterations: break
      i+=1
      print("nans",np.sum(nans))
      boundary=find_boundaries(~nans, mode='outer')
      tmp = neighborhood_average(ret, neighborhood)
      ret[boundary] = tmp[boundary]
      ret[mask] = grid[mask]
      nans = np.isnan(ret)
    return ret

gaussian_kernel(voxels=5, sig=1.0, ndim=3)

creates gaussian kernel with side length l and a sigma of sig

Source code in mrdna/arbdmodel/submodule/grid.py
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
def gaussian_kernel(voxels=5, sig=1., ndim=3):
    ## TODO: rewrite using isotropic_kernel
    """
    creates gaussian kernel with side length `l` and a sigma of `sig`
    """
    gauss_list = []
    try:
        sig[0]
    except:
        sig = [sig]*ndim
    try:
        voxels[0]
    except:
        voxels = [voxels]*ndim

    for l,s in zip(voxels,sig):
        ax = np.linspace(-(l - 1) / 2., (l - 1) / 2., l)
        gauss_list.append( np.exp(-0.5 * np.square(ax) / np.square(s)) )

    gauss_list2 = []
    kernel = gauss_list[0] # np.empty((0))
    for i,g in enumerate(gauss_list[1:]):
        i=i+1
        sl = tuple(slice(None) if j==i else None for j in range(i))
        kernel = kernel[...,None]*g[sl]

    return kernel / np.sum(kernel)

isotropic_kernel(function, delta, shape=None, cutoff=None, normalize=False)

Generate an isotropic kernel function.

Parameters:

Name Type Description Default
function callable

A callable representing the kernel function. It should accept a scalar or an array-like object and return a scalar or an array-like object of the same shape.

required
delta list or tuple

A list or tuple specifying the spacing between grid points in each dimension.

required
shape list or tuple

A list or tuple specifying the shape of the kernel. If None, it is determined based on the cutoff parameter.

None
cutoff scalar

A scalar specifying the maximum distance from the center to consider for the kernel. If provided, it determines the shape of the kernel based on the delta values.

None
normalize bool

A boolean indicating whether to normalize the kernel. If True, the kernel values are divided by the sum of all kernel values, ensuring that the kernel sums to 1.

False

Returns:

Name Type Description
ndarray

An ndarray representing the isotropic kernel.

Raises:

Type Description
ValueError

If the cutoff parameter is None, as it is required to determine the shape of the kernel.

Example

import numpy as np def gaussian(x): ... return np.exp(-0.5 * x**2) delta = [0.1, 0.1, 0.1] kernel = isotropic_kernel(gaussian, delta, cutoff=1.0, normalize=True) print(kernel)

Source code in mrdna/arbdmodel/submodule/grid.py
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
def isotropic_kernel( function, delta, shape = None, cutoff = None, normalize = False):
    """
    Generate an isotropic kernel function.

    Parameters:
        function (callable): A callable representing the kernel function.
            It should accept a scalar or an array-like object and return a scalar or an array-like object of the same shape.
        delta (list or tuple): A list or tuple specifying the spacing between grid points in each dimension.
        shape (list or tuple, optional): A list or tuple specifying the shape of the kernel. If None, it is determined based on the cutoff parameter.
        cutoff (scalar, optional): A scalar specifying the maximum distance from the center to consider for the kernel. If provided, it determines the shape of the kernel based on the delta values.
        normalize (bool, optional): A boolean indicating whether to normalize the kernel. If True, the kernel values are divided by the sum of all kernel values, ensuring that the kernel sums to 1.

    Returns:
        ndarray: An ndarray representing the isotropic kernel.

    Raises:
        ValueError: If the cutoff parameter is None, as it is required to determine the shape of the kernel.

    Example:
        >>> import numpy as np
        >>> def gaussian(x):
        ...     return np.exp(-0.5 * x**2)
        >>> delta = [0.1, 0.1, 0.1]
        >>> kernel = isotropic_kernel(gaussian, delta, cutoff=1.0, normalize=True)
        >>> print(kernel)

    """

    if cutoff is not None:
        shape = [(cutoff*2)//d for d in delta]
    elif cutoff is None:
        raise ValueError

    ndim = len(shape)
    assert( len(delta) == ndim )

    centers = [(np.arange(n)-(n-1)*0.5)*d for n,d in zip(shape,delta)]
    CENTERS = np.meshgrid(*centers, indexing='ij')
    R = np.sqrt( sum([C**2 for C in CENTERS]) )
    kernel = function(R)
    if normalize:
        kernel = kernel/np.sum(kernel)
    return kernel

neighborhood_average(grid, neighborhood=1, fill_value='mirror')

Compute the neighborhood average of a grid.

Parameters:

Name Type Description Default
grid ndarray

The input grid.

required
neighborhood int or list

The size of the neighborhood to use for averaging. If an integer is provided, the same neighborhood size is used for all dimensions. If a list is provided, it should specify the neighborhood size for each dimension separately. Default is 1.

1
fill_value str or ndarray

The fill value to use for padding the grid. If 'mirror', values are mirrored along each dimension. If 'nearest', values are filled with the nearest neighbor. If 'periodic', values are filled with periodic boundary conditions. If a number is provided, it is used directly as the fill value. Default is 'mirror'.

'mirror'

Returns:

Name Type Description
ndarray

The grid with the neighborhood average computed.

Notes
  • Currently, only the 'mirror' fill value option is implemented.
  • This function requires the 'average_grids' function from an external source.
Example

grid = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) averaged_grid = neighborhood_average(grid, neighborhood=2, fill_value='mirror') print(averaged_grid)

Source code in mrdna/arbdmodel/submodule/grid.py
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def neighborhood_average(grid, neighborhood = 1, fill_value='mirror'):
    """
    Compute the neighborhood average of a grid.

    Parameters:
        grid (ndarray): The input grid.
        neighborhood (int or list): The size of the neighborhood to use for averaging. If an integer is provided, the same neighborhood size is used for all dimensions. If a list is provided, it should specify the neighborhood size for each dimension separately. Default is 1.
        fill_value (str or ndarray): The fill value to use for padding the grid. If 'mirror', values are mirrored along each dimension. If 'nearest', values are filled with the nearest neighbor. If 'periodic', values are filled with periodic boundary conditions. If a number is provided, it is used directly as the fill value. Default is 'mirror'.

    Returns:
        ndarray: The grid with the neighborhood average computed.

    Notes:
        - Currently, only the 'mirror' fill value option is implemented.
        - This function requires the 'average_grids' function from an external source.

    Example:
        >>> grid = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
        >>> averaged_grid = neighborhood_average(grid, neighborhood=2, fill_value='mirror')
        >>> print(averaged_grid)
    """
    try:
        neighborhood[0]
    except:
        neighborhood = [neighborhood] * grid.ndim

    padded_grid = np.ones( [s+dx*2 for s,dx in zip(grid.shape,neighborhood)] )*np.nan
    if fill_value == 'mirror':
      def get_slices(mirror_dims, left):
        sl1 = tuple(slice(n,-n) if i not in mirror_dims else 
                    slice(None,n) if left[mirror_dims.index(i)]
                    else slice(-n,None)
                    for i,n in enumerate(neighborhood))
        sl2 = tuple(slice(None,None) if i not in mirror_dims else
                    slice((n-1),None,-1) if left[mirror_dims.index(i)]
                    else slice(None,-(n+1),-1)
                    for i,n in enumerate(neighborhood))
        return sl1,sl2

      def build_dim_arrays(nd, mirror_dims, left):
        d0 = mirror_dims[-1] if len(mirror_dims) > 0 else -1
        mirror_dims.append(d0)
        left.append(False)
        for d in range(d0+1,grid.ndim+1-nd):
          mirror_dims[-1] = d
          left[-1] = False
          yield list(mirror_dims), list(left)
          left[-1] = True
          yield list(mirror_dims), list(left)

      all_mirror_dims_arrays = []
      all_left_arrays = []
      for nd in range(1,grid.ndim+1):
        # print(nd,grid.ndim)
        mirror_dims_arrays = [[]]
        left_arrays = [[]]
        while nd > 0:
          for m0,l0 in list(zip(mirror_dims_arrays,left_arrays)):
            for m,l in build_dim_arrays(nd, m0, l0):
              mirror_dims_arrays.append(m)
              left_arrays.append(l)
          nd=nd-1
        # print("dbg",nd)  
        # print(mirror_dims_arrays)
        all_mirror_dims_arrays.extend(mirror_dims_arrays)
        all_left_arrays.extend(left_arrays)

      for m,l in zip(all_mirror_dims_arrays, all_left_arrays):
        ## set values mirroring d
        sl1,sl2 = get_slices(m, l)
        padded_grid[sl1] = grid[sl2]
    elif fill_value == 'nearest':
      raise NotImplementedError      
    elif fill_value == 'periodic':
      raise NotImplementedError
    else:
      padded_grid = fill_value   

    sl = tuple(slice(n,-n) for n in neighborhood)
    padded_grid[sl] = grid
    # print(padded_grid)

    ## Gather sliced arrays to perform average
    def get_slice(neighborhood, slices=None):
      # print("get_slice",neighborhood)
      if slices is None: slices = [[]]
      if len(neighborhood) > 1:
        slices = get_slice(neighborhood[1:], slices)
      nmax = 2*neighborhood[0]
      slices = [[slice(n,-(nmax-n) if nmax != n else None)]+sls for n in range(nmax+1) for sls in slices]
      return slices

    slices = get_slice(neighborhood)
    result = average_grids([padded_grid[sl] for sl in slices])
    return result