CODEWITHSUNDEEP
pythonnumpyscipyscikit-learnastropy
Dalek
Dalek

Reputation: 4318

Determine the coordinates of local maximas in a two-dimensional array using derivative

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.

Upvotes: 1

Views: 1177

Answers (2)

Christoph
Christoph

Reputation: 2940

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.

Upvotes: 3

kezzos
kezzos

Reputation: 3221

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.

Upvotes: 4

Related Questions