1

I have a fits image and I am trying to find the coordinates of local maxima in my image but so far I couldn't quite make it work. My image can be find here. What I have so far is 

import numpy as np
import scipy.nimage as ndimage
from astropy.wcs import WCS
from astropy import units as u
from astropy import coordinates as coord
from astropy.io import fits
import scipy.ndimage.filters as filters
from scipy.ndimage.filters import maximum_filter
hdulist=fits.open("MapSNR.fits")
#reading a two dimensional array from fits file
d=hdulist[0].data
w=WCS("MapSNR.fits") 
idx,idy=np.where(d==np.max(d))
rr,dd=w.all_pix2word(idx,idy,o)
c=coord.SkyCoord(ra=rr*u.degree, dec=dd*u.degree)
#The sky coordinate of the image maximum
print c.ra
print c.dec

That is how I can find the global maximum of the image, but I would like to obtain the coordinates of the local maximas which have the significance of being greater than three.

What I have found by looking up in the web was this following answer which doesn't work properly in my case. update: I have used this function

def detect_peaks(data, threshold=1.5, neighborhood_size=5):

  data_max = filters.maximum_filter(data, neighborhood_size)
  maxima = (data == data_max)
  data_min = filters.minimum_filter(data, neighborhood_size)
  diff = ((data_max - data_min) > threshold)
  maxima[diff == 0] = 0 # sets values <= threshold as background
  labeled, num_objects = ndimage.label(maxima)
  slices = ndimage.find_objects(labeled)
  x,y=[],[]
  for dy,dx in slices:
    x_center = (dx.start + dx.stop - 1)/2
    y_center = (dy.start + dy.stop - 1)/2
    x.append(x_center)
    y.append(y_center)
  return x,y

I would like to find a method using a better approach like the derivative in the array or divide and conquer method. I will appropriate for a better recommended solution.

Community
  • 1
  • 1
Dalek
  • 4,168
  • 11
  • 48
  • 100
  • Why does the method you describe (the answer) fail in your case? – kezzos Aug 16 '15 at 20:39
  • @kezzos it doesn't return the right values.. – Dalek Aug 16 '15 at 20:43
  • Could you explain why the answer at http://stackoverflow.com/questions/9111711/get-coordinates-of-local-maxima-in-2d-array-above-certain-value "doesn't work properly"? – Warren Weckesser Aug 16 '15 at 23:49
  • @WarrenWeckesser I have tried different threshold values but at the end the indices the code returns do not corresond to the right maximas. – Dalek Aug 16 '15 at 23:55
  • Are you using a simple definition of a local maximum? I.e. a point where the value is greater than or equal to the values of its eight neighbors? – Warren Weckesser Aug 17 '15 at 00:02
  • ... or are you trying to avoid lots of "uninteresting" local maxima that can arise in a noisy image? – Warren Weckesser Aug 17 '15 at 00:04
  • @WarrenWeckesser I would like to obtain local maximas with values greater than 3, meaning I also do not want to obtain neighboring points with more and less the same values. – Dalek Aug 17 '15 at 00:08
  • possible duplicate of [Peak detection in a 2D array](http://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array) – Iguananaut Aug 17 '15 at 19:57

2 Answers2

4

So I have this, using skimage adaptive threshold. Hope it helps:

Original enter image description here

Code

from skimage.filters import threshold_adaptive
import matplotlib.pyplot as plt
from scipy import misc, ndimage
import numpy as np

im = misc.imread('\Desktop\MapSNR.jpg')

# Apply a threshold
binary_adaptive = threshold_adaptive(im, block_size=40, offset=-20).astype(np.int)
# Label regions and find center of mass
lbl = ndimage.label(binary_adaptive)
points = ndimage.measurements.center_of_mass(binary_adaptive, lbl[0], [i+1 for i in range(lbl[1])])

for i in points:
    p = [int(j) for j in i]
    binary_adaptive[i] += 5

plt.figure()
plt.imshow(im, interpolation='nearest', cmap='gray')
plt.show()

plt.figure()
plt.imshow(binary_adaptive, interpolation='nearest', cmap='gray')
plt.show()

Output

enter image description here

Changing the parameters of the threshold will have a large effect on where the local maxima will be found, and how many maxima are found.

kezzos
  • 3,023
  • 3
  • 20
  • 37
3

You could use the photutils.detection.find_peaks function, which is one of the photutils detection methods.

If you look at the photutils.detection.find_peaks implementation, you'll see that it's using scipy.ndimage.maximum_filter to compute a maximum image (by default in a 3x3 box size footprint) and finding the pixels where the original image equals the maximum image.

The rest of the function is mostly for two things that might also be of interest to you:

  1. if you pass in a wcs object, you can get sky coordinates out, not just pixel coordinates.
  2. there's an option to get sub-pixel precision coordinates.
Christoph
  • 2,790
  • 2
  • 18
  • 23