title: Peak detection in a noisy 2d array
body: <p>I'm trying to get <code>python</code> to return, as close as possible, the <strong>center</strong> of the most obvious clustering in an image like the one below:</p>
  
  <p><img src="https://i.stack.imgur.com/iE0C3.png" alt="image"></p>
  
  <p>In my <a href="https://stackoverflow.com/questions/16822334/find-peak-of-2d-histogram">previous question</a> I asked how to get the global maximum and the local maximums of a 2d array, and the answers given worked perfectly. The issue is that the center estimation I can get by averaging the global maximum obtained with different bin sizes is always slightly off than the one I would set <em>by eye</em>, because I'm only accounting for the biggest <strong>bin</strong> instead of a <strong>group</strong> of biggest bins (like one does by eye).</p>
  
  <p>I tried adapting the <a href="https://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array">answer to this question</a> to my problem, but it turns out my image is <em>too noisy</em> for that algorithm to work. Here's my code implementing that answer:</p>
  
  <pre><code>import numpy as np
  from scipy.ndimage.filters import maximum_filter
  from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
  import matplotlib.pyplot as pp
  
  from os import getcwd
  from os.path import join, realpath, dirname
  
  # Save path to dir where this code exists.
  mypath = realpath(join(getcwd(), dirname(__file__)))
  myfile = 'data_file.dat'
  
  x, y = np.loadtxt(join(mypath,myfile), usecols=(1, 2), unpack=True)
  xmin, xmax = min(x), max(x)
  ymin, ymax = min(y), max(y)
  
  rang = [[xmin, xmax], [ymin, ymax]]
  paws = []
  
  for d_b in range(25, 110, 25):
      # Number of bins in x,y given the bin width 'd_b'
      binsxy = [int((xmax - xmin) / d_b), int((ymax - ymin) / d_b)]
  
      H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
      paws.append(H)
  
  
  def detect_peaks(image):
      """
      Takes an image and detect the peaks usingthe local maximum filter.
      Returns a boolean mask of the peaks (i.e. 1 when
      the pixel's value is the neighborhood maximum, 0 otherwise)
      """
  
      # define an 8-connected neighborhood
      neighborhood = generate_binary_structure(2,2)
  
      #apply the local maximum filter; all pixel of maximal value 
      #in their neighborhood are set to 1
      local_max = maximum_filter(image, footprint=neighborhood)==image
      #local_max is a mask that contains the peaks we are 
      #looking for, but also the background.
      #In order to isolate the peaks we must remove the background from the mask.
  
      #we create the mask of the background
      background = (image==0)
  
      #a little technicality: we must erode the background in order to 
      #successfully subtract it form local_max, otherwise a line will 
      #appear along the background border (artifact of the local maximum filter)
      eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)
  
      #we obtain the final mask, containing only peaks, 
      #by removing the background from the local_max mask
      detected_peaks = local_max - eroded_background
  
      return detected_peaks
  
  
  #applying the detection and plotting results
  for i, paw in enumerate(paws):
      detected_peaks = detect_peaks(paw)
      pp.subplot(4,2,(2*i+1))
      pp.imshow(paw)
      pp.subplot(4,2,(2*i+2) )
      pp.imshow(detected_peaks)
  
  pp.show()
  </code></pre>
  
  <p>and here's the result of that (varying the bin size):</p>
  
  <p><img src="https://i.stack.imgur.com/duZ7f.png" alt="enter image description here"></p>
  
  <p>Clearly my background is too noisy for that algorithm to work, so the question is: <strong>how can I make that algorithm less sensitive</strong>? If an alternative solution exists then please let me know.</p>
  
  <hr>
  
  <h2>EDIT</h2>
  
  <p>Following Bi Rico advise I attempted smoothing my 2d array before passing it on to the local maximum finder, like so:</p>
  
  <pre><code>H, xedges, yedges = np.histogram2d(x, y, range=rang, bins=binsxy)
  H1 = gaussian_filter(H, 2, mode='nearest')
  paws.append(H1)
  </code></pre>
  
  <p>These were the results with a <code>sigma</code> of 2, 4 and 8:</p>
  
  <p><img src="https://i.stack.imgur.com/ZciZY.png" alt="enter image description here"></p>
  
  <h2>EDIT 2</h2>
  
  <p>A <code>mode ='constant'</code> seems to work much better than <code>nearest</code>. It converges to the right center with a <code>sigma=2</code> for the largest bin size:</p>
  
  <p><img src="https://i.stack.imgur.com/AYc24.png" alt="enter image description here"></p>
  
  <p>So, <strong>how do I get the coordinates of the maximum that shows in the last image?</strong></p>