Efficient density function computation

Question

I have a large image in numpy array form (opencv returns it as a 2d array of 3 uint8 values) and want to compute a sum of gaussian kernels for each pixel, i.e. (there's still no LaTeX support in SO is there?):

for N different kernels with a specified weight w, mean and diagonal covariance matrix.

So basically I want a function compute_densities(image, kernels) -> numpy array of floats. What's the best way to do this efficiently in python? I'd be surprised if there wasn't already a library function in scipy for this, but I had statistics at uni a long time ago, so I do get a bit confused with the details of the documentation..

Basically I want the following, just way more efficient than naive python (2pi^{-3/2} is ignored since it's a constant factor that doesn't matter for me since I'm only interested in ratios between the probabilities)

def compute_probabilities(img, kernels):
    np.seterr(divide='ignore') # 1 / covariance logs an error otherwise
    result = np.zeros((img.shape[0], img.shape[1]))
    for row_pos, row_val in enumerate(img):
        for col_pos, val in enumerate(row_val):
            prob = 0.0
            for kernel in kernels:
                mean, covariance, weight = kernel
                val_sub_mu = np.array([val]).T - mean
                cov_inv = np.where(covariance != 0, 1 / covariance, 0)
                tmp = val_sub_mu.T.dot(cov_inv).dot(val_sub_mu)
                prob += weight / np.sqrt(np.linalg.norm(covariance)) * \
                        math.exp(-0.5 * tmp)
            result[row_pos][col_pos] = prob
    np.seterr(divide='warn')
    return result

Input: cv2.imread on some jpg, which gives a 2d array (height x width) of a 3 uint8 struct containing the 3 color channels.

Kernels is a namedtuple('Kernel', 'mean covariance weight'), mean is a vector, covariance is a 3x3 matrix with everything but the diagonal being zero and weight is a float 0 < weight < 1. For simplicity I only specify the diagonals and then convert it to a 3x3 matrix afterwards: (the representation isn't set in stone I don't care how it's represented so be free to change all of that):

some_kernels = [
   Kernel(np.array([(73.53, 29.94, 17.76)]), np.array([(765.40, 121.44, 112.80)]), 0.0294),
   ...
]

def fixup_kernels(kernels):
    new_kernels = []
    for kernel in kernels:
        cov = np.zeros((3, 3))
        for pos, c in enumerate(kernel.covariance[0]):
            cov[pos][pos] = c
        new_kernels.append(Kernel(kernel.mean.T, cov, kernel.weight))
    return new_kernels

 some_kernels = fixup_kernels(some_kernels)
 img = cv2.imread("something.jpg")
 result = compute_probabalities(img, some_kernels)

Answer 1

EDIT

I verified that this produces same results as the original code:

def compute_probabilities_fast(img, kernels):
    np.seterr(divide='ignore')
    result = np.zeros((img.shape[0], img.shape[1]))
    for kernel in kernels:
        mean, covariance, weight = kernel
        cov_inv = np.where(covariance != 0, 1 / covariance, 0)
        mean = mean[:,0]
        img_sub_mu = img - mean
        img_tmp = np.sum( img_sub_mu.dot(cov_inv) * img_sub_mu, axis=2 )
        result += (weight / np.sqrt(np.linalg.norm(covariance))) * np.exp(-0.5 * img_tmp)
    return result

Explanation:

mean[:,0] makes the shape simply (3,) instead of (3,1).

img - mean broadcasts to the whole image and subtracts the means from each pixel.

img_sub_mu.dot(cov_inv) is roughly equivalent to val_sub_mu.T.dot(cov_inv).

np.sum( ... * img_sub_mu, axis=2 ) is roughly equivalent to .dot(val_sub_mu). Can't use dot, though, because doing that will add extra dimensions. For example an array M x N x K dotted with an array M x K x N would produce result M x N x M x N, dot behaves differently on one-dimensional and multi-dimensional data. So we just do an element-wise multiply and then sum along the last dimension.

Actually the "sum of gaussian kernels" bit in the question had me confused. What is requested is a calculation in which, for each output pixel the value depends only on the input value for the same pixel, but not on the values of neighboring pixels. So, this is nothing like a gaussian blur (which would use convolution), it is just a calculation performed on each pixel individually.

P.S. 1 / covariance is problematic. Are you sure you don't want np.linalg.inv(covariance) instead?

OLD ANSWER

It sounds like what you want is one of these:

scipy.signal.convolve2d

scipy.ndimage.filters.convolve

The question is a bit confusing, are you trying to calculate a bunch of images convolved with different gaussians, or a single image convolved with a sum of gaussians? Is your kernel separable? (If yes, use two convolutions Mx1 and 1xN instead of one MxN) The scipy functions you'd use are the same in any case.

Also of course you'd want to pre-calculate your kernels with a combination of numpy.random.normal and meshgrid.

Answer 2

As I understand it, your goal is to evaluate the probability density of each image pixels value with respect to a mixture of multivariate normal distributions.

[There is currently (2013-11-20) a mistake in the code in your question and @Alex I's answer - the | | surrounding \Sigma in the above equation actually denotes the determinant rather than the vector norm - see e.g. here. In the case of a diagonal covariance the determinant is just the product of the diagonal elements.]

It is possible to implement the density calculation very efficiently in terms of numpy array operations. The following implementation makes use of the spherical (i.e. diagonal) nature of the covariance matrices in your problem:

def compute_probabilities_faster(img, kernels):
  means, covs, weights = map(np.dstack, zip(*kernels)) 
  pixels_as_rows = img.reshape((-1, 3, 1))
  responses = np.exp(-0.5 * ((pixels_as_rows - means) ** 2 / covs).sum(axis=1))
  factors = 1. / np.sqrt(covs.prod(axis=1) * ((2 * np.pi) ** 3))
  return np.sum(responses * factors * weights, axis=2).reshape(img.shape[:2])

This function operates directly on the kernels as they are initially represented, ie. without modification by your fixup_kernels function. When the normalisation factor (2 * np.pi) ** 3 is removed (and calls to linalg.norm are replaced with linalg.det) this function matches the output of your code (well enough to satisfy np.allclose).

The closest out-of-the-box functionality in SciPy (as of 0.13) is the implementation of Kernel Density Estimation in scipy.stats (see here) which defines a very similar distribution where the covariance matrices are identical for each kernel - for this reason not appropriate for your problem.