Databases
StackOverflow2013
Postsdbo
POPeak detection in a noisy 2d array
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POPeak detection in a noisy 2d array
    primarykey
    data
    text
    <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>
    singulars
    1. PO
      singulars
      1. PTAnswer
    1. This table or related slice is empty.
    1. PTQuestion
    1. USCommunity
    1. USGabriel
    plurals
    1. PL
      singulars
      1. LTLinked
    2. PL
      singulars
      1. LTLinked
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    2. PO
      singulars
      1. PTAnswer
    3. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. POPeak detection in a noisy 2d array
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. POPeak detection in a noisy 2d array
      1. This table or related slice is empty.
      1. VTUpMod
    3. VO
      singulars
      1. POPeak detection in a noisy 2d array
      1. USnicholaschris
      1. VTFavorite
    1. COHave you tried smoothing your data before applying your algorithm? A Gaussian and/or median filters might help.
      singulars
      1. POPeak detection in a noisy 2d array
      1. USBi Rico
    2. COSee updated question please.
      singulars
      1. POPeak detection in a noisy 2d array
      1. USGabriel
    3. COHow about, `np.unravel_index(array.argmax(), array.shape)`.
      singulars
      1. POPeak detection in a noisy 2d array
      1. USBi Rico
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload