Source code for pysit.util.image_processing

import itertools

import numpy as np
import scipy
import scipy.signal as signal

__all__ = ['gaussian_kernel', 'blur_image', 'resample_array']

[docs]def gaussian_kernel(size, sigma, mesh_deltas=None): """ Returns a normalized gauss kernel array for convolutions. Parameters ---------- size : iterable Number of pixels per standard deviation sigma : float or iterable Standard deviation of gaussian in physical units mesh_deltas : iterable, optional Grid size in physical units. Defaults to array of ones (pixel units). """ imsize = size.astype(int) if mesh_deltas is None: mesh_deltas = np.ones(imsize) # Use numpy index tricks to generate a meshgrid like array for each dimension. # Product with the delta ensures that the indices are in the physical units # specified. idx = tuple([slice(-s, s+1) for s in imsize]) grid = [g.astype(np.float)*dx for g,dx in zip(np.mgrid[idx], mesh_deltas)] # Create a different kernel if sigma is specifed differently for different axes. if np.iterable(sigma) and len(sigma) == len(grid): arg = sum([dim.astype(float)**2/(2*(s**2)) for dim,s in zip(grid,sigma)]) else: arg = sum([dim**2/(2*float(sigma)**2) for dim in grid]) g = np.exp(-(arg)) # Return the normalized Gaussian return g/g.sum()
[docs]def blur_image(im, sigma=None, freq=None, mesh_deltas=None, n_sigma=1.0): """ Returns a blurred image by convolving with a gaussian kernel. Parameters ---------- im : ndarray Input image sigma : float or ndarray Standard deviation of blur kernel freq : Cutoff (spatial) frequency of the blur kernel. Equivalent to 1/sigma. mesh_deltas : iterable, optional Grid size in physical units. Defaults to array of ones (pixel units). n_sigma : float Width of the blur kernel in standard deviations. Notes ----- * One of `sigma` or `freq` must be specified. If both are specified, `freq` takes priority. """ if sigma is None and freq is None: raise ValueError('Either sigma or frequency must be set.') if mesh_deltas is None: mesh_deltas = np.ones(im.ndim) if not np.iterable(mesh_deltas): mesh_deltas = np.array([mesh_deltas]) else: mesh_deltas = np.asarray(mesh_deltas) # frequency takes priority if both freq and sigma are specified. if freq is not None: sigma = 1./freq # determine the size, in pixels, of the kernel kernel_size_pixel = (np.ceil(n_sigma*sigma / mesh_deltas)).astype('int32') kernel = gaussian_kernel(kernel_size_pixel, sigma, mesh_deltas) # Pad the image by the kernel size on all sides, to get smoothing without # edge effects im_ = np.pad(im, list(zip(kernel_size_pixel,kernel_size_pixel)), mode='edge') improc = signal.fftconvolve(im_, kernel, mode='valid') #return a copy to ensure a contiguous array return improc.copy()
[docs]def resample_array(arr, new_size, mode='nearest'): """ Returns a resampled array at new resolution. Parameters ---------- arr : ndarray Input array new_size : iterable Size of the output array mode : {'nearest', 'linear'} Interpolation mode. """ sh = arr.shape y = arr # For each dimension, construct a 1D interpolator and interpolate. # new size, old size for nsz, osz, i in zip(new_size, sh, itertools.count()): x = np.arange(osz, dtype=np.float) I = scipy.interpolate.interp1d(x, y, copy=False, axis=i, kind=mode) new_x = np.linspace(0, osz-1, nsz) y = I(new_x) return y