find_peaks.py

# shamelessly hacked for Albert Cardona's example

from ij import WindowManager
from ij.plugin.frame import RoiManager
from ij.gui import PointRoi
from script.imglib.analysis import DoGPeaks
from script.imglib.color import Red
from script.imglib.algorithm import Scale3D
from script.imglib.math import Compute
from script.imglib import ImgLib
from ij3d import Image3DUniverse
from javax.vecmath import Color3f, Point3f
from ij import IJ

cell_diameter = 4.0 # in microns
minPeak = 20.0 # The minimum intensity for a peak to be considered so.


imp = IJ.getImage()

# Scale the X,Y axis down to isotropy with the Z axis
cal = imp.getCalibration()
scale2D = cal.pixelWidth / cal.pixelDepth

print scale2D
#iso = Compute.inFloats(Scale2D(ImgLib.wrap(imp), scale2D))
# using scale3D shold allow me to scale z in the future, instead of only downscaling XY
iso = Compute.inFloats(Scale3D(ImgLib.wrap(imp), scale2D, scale2D,1 ))
print type(iso)


# Find peaks by difference of Gaussian
sigma = (cell_diameter  / cal.pixelWidth) * scale2D
peaks = DoGPeaks(iso, sigma, sigma * 0.5, minPeak, 1).getPeaks()
print "Found", len(peaks), "peaks"

# Only keep maxima (which are minima in DoG)
maxima = []
for p in peaks:
  if p.isMin():
    maxima.append(p)

# Convert the peaks into points in calibrated image space
ps = []
if "ROI Manager" in WindowManager.getNonImageTitles():
    WindowManager.getWindow("ROI Manager").close()

rm = RoiManager()
rm.show()


for peak in maxima:
  cellRoi = PointRoi(peak.getPosition(0)/scale2D,peak.getPosition(1)/scale2D)
  cellRoi.setPosition(int(peak.getPosition(2)))
  rm.addRoi(cellRoi)

print "Found", len(maxima), "peaks"

print type(iso)
# ImgLib.wrap(iso).show()

# Show the peaks as spheres in 3D, along with orthoslices:
#univ = Image3DUniverse(1024, 1024)
#univ.addIcospheres(ps, Color3f(1, 0, 0), 2, cell_diameter/2, "Cells").setLocked(True)
#univ.addOrthoslice(imp).setLocked(True)
#univ.show()