title: How to efficiently compute the heat map of two Gaussian distribution in Python?
body: <p>I am trying to produce a heat map where the pixel values are governed by two independent 2D Gaussian distributions. Let them be Kernel1 (muX1, muY1, sigmaX1, sigmaY1) and Kernel2 (muX2, muY2, sigmaX2, sigmaY2) respectively. To be more specific, the length of each kernel is three times its standard deviation. The first Kernel has sigmaX1 = sigmaY1 and the second Kernel has sigmaX2 &lt; sigmaY2. Covariance matrix of both kernels are diagonal (X and Y are independent). Kernel1 is usually completely inside Kernel2. </p>
  
  <p>I tried the following two approaches and the results are both unsatisfactory. Can someone give me some advice?</p>
  
  <p>Approach1:</p>
  
  <p>Iterate over all pixel value pairs (i, j) on the map and compute the value of I(i,j) given by I(i,j) = P(i, j | Kernel1, Kernel2) = 1 - (1 - P(i, j | Kernel1)) * (1 - P(i, j | Kernel2)). Then I got the following result, which is good in terms of smoothness. But it takes 10 seconds to run on my computer, which is too slow. </p>
  
  <p>Codes:</p>
  
  <pre><code>def genDensityBox(self, height, width, muY1, muX1, muY2, muX2, sigmaK1, sigmaY2, sigmaX2):
      densityBox = np.zeros((height, width))
      for y in range(height):
          for x in range(width):
              densityBox[y, x] += 1. - (1. - multivariateNormal(y, x, muY1, muX1, sigmaK1, sigmaK1)) * (1. - multivariateNormal(y, x, muY2, muX2, sigmaY2, sigmaX2))
      return densityBox
  
  def multivariateNormal(y, x, muY, muX, sigmaY, sigmaX):
      return norm.pdf(y, loc=muY, scale=sigmaY) * norm.pdf(x, loc=muX, scale=sigmaX)
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/l9g3m.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/l9g3m.png" alt="First approach"></a></p>
  
  <p>Approach2:</p>
  
  <p>Generate two images corresponding to two kernels separately and then blend them together using certain alpha value. Each image is generated by taking the outer product of two one-dimensional Gaussian filter. Then I got the following result, which is very crude. But the advantage of this approach is that it is very fast due to the use of outer product between two vectors.
  <a href="https://i.stack.imgur.com/riaMq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/riaMq.png" alt="Second approach"></a></p>
  
  <p>Since the first one is slow and the second on is crude, I am trying to find a new approach that achieves good smoothness and low time-complexity at the same time. Can someone give me some help?</p>
  
  <p>Thanks!</p>
  
  <p>For the second approach, the 2D Gaussian map can be easily generated as mentioned <a href="https://stackoverflow.com/questions/29731726/how-to-calculate-a-gaussian-kernel-matrix-efficiently-in-numpy">here</a>:</p>
  
  <pre><code>def gkern(self, sigmaY, sigmaX, yKernelLen, xKernelLen, nsigma=3):
      """Returns a 2D Gaussian kernel array."""
      yInterval = (2*nsigma+1.)/(yKernelLen)
      yRow = np.linspace(-nsigma-yInterval/2.,nsigma+yInterval/2.,yKernelLen + 1)
      kernelY = np.diff(st.norm.cdf(yRow, 0, sigmaY))
      xInterval = (2*nsigma+1.)/(xKernelLen)
      xRow = np.linspace(-nsigma-xInterval/2.,nsigma+xInterval/2.,xKernelLen + 1)
      kernelX = np.diff(st.norm.cdf(xRow, 0, sigmaX))    
      kernelRaw = np.sqrt(np.outer(kernelY, kernelX))
      kernel = kernelRaw / (kernelRaw.sum())
      return kernel
  </code></pre>