Posted by
Michael Schmid on
Mar 19, 2007; 6:55pm
URL: http://imagej.273.s1.nabble.com/ImageJ-Findpeak-plugin-tp3700022p3700025.html
Hi Josh,
unfortunately the plugin works with almost perfect periodicity
only. It does not search for the maxima of the line itself but
rather for the maxima of the autocorrelation along the line.
You can make it less sensitive to varying periodicities by
removing the line "if (Double.isNaN(x)) return 0;" and by
setting it to search for the shortest, not for the best
periodicity.
It won't give you any statistics about the distribution of
maxima.
It is pretty long (24 kB); I hope that the mailinglist won't
truncate it.
Michael
________________________________________________________________
On 19 Mar 2007, at 16:02, spinodalman wrote:
> Yes, my final aim is to mesaure the periodicity of an array of
> features.
> These features are not ususally perfectly periodic, however, and I
> would
> like to get some statistical information about the periodicity as
> well. I am
> very interested in the plugin.
>
> Josh
>
________________________________________________________________
import ij.*;
import ij.plugin.frame.PlugInFrame;
import ij.process.*;
import ij.measure.*;
import ij.gui.*;
import ij.util.Tools;
import java.awt.*;
import java.awt.event.*;
import java.util.*;
/**
This plugin continuously measures the periodicity along a line scan
*/
public class Measure_Periodicity extends PlugInFrame
implements MouseListener, MouseMotionListener, KeyListener,
ActionListener, ItemListener, Runnable {
//MouseListener, MouseMotionListener, KeyListener to watch an
ImagePlus
//KeyListener for text field (go to...)
//ActionListener for Buttons and textfields (number input)
//ItemListener for Choice (Pulldown Menus)
//Runnable for background thread
private static Measure_Periodicity instance; //we want only
one instance of this; remember it
private ImagePlus imp, lastImp; //the ImagePlus
that we listen to and the last one
private ImagePlus plotImage;
private Thread bgThread; //thread for
plotting (in the background)
private final static int PROFILE = 1;
private final static int CORRELATION = 2;
private int showPlot = 0; //could be 0 (no
plot), PROFILE, CORRELATION
private double angle;
private double periodicity;
private double result;
private String xUnit;
private double xInc;
double[] yProfile;
double[] yCorr;
private Button updateB;
private Button toResultB;
private TextField[] numberT;
private Choice methodChoice;
private Choice plotChoice;
private Choice[] functionChoice;
private Label imageLabel;
private Label statusLabel;
private Label periodLabel;
private Label resultLabel;
/** The constructor, preparing the window (panel) */
public Measure_Periodicity() {
super("Measure Periodicity");
if (instance!=null) { //we want only
one instance
instance.toFront();
return;
}
instance = this; //we are new,
set up the window
WindowManager.addWindow(this);
imageLabel = new Label( //set up all
the panel items
"No Image "); //labels with
enough width for later contents
statusLabel = new Label("No Line Selection");
periodLabel = new Label("<no value yet>");
resultLabel = new Label("<no value yet>");
int trim = IJ.isMacOSX()?11:0;
updateB = new TrimmedButton("Image", trim);
toResultB = new TrimmedButton("to Results", trim);
numberT = new TextField[2];
numberT[0] = new TextField("1.00",10);
numberT[1] = new TextField("1.00",10);
String[] methods = new String[] {"Best
Periodicity","Shortest Periodicity"};
String[] plots = new String[] {"Plot Nothing", "Plot
Profile", "Plot Autocorrelation"};
String[] functions = new String[] {"*","/","* sqrt", "/sqrt"};
methodChoice = new Choice();
plotChoice = new Choice();
functionChoice = new Choice[2];
functionChoice[0] = new Choice();
functionChoice[1] = new Choice();
for (int i=0; i<methods.length; i++)
methodChoice.add(methods[i]);
for (int i=0; i<plots.length; i++)
plotChoice.add(plots[i]);
for (int i=0; i<functions.length; i++)
for (int j=0; j<functionChoice.length; j++)
functionChoice[j].add(functions[i]);
updateB.addActionListener(this); //set up the
listeners (callback routines) for input
toResultB.addActionListener(this);
numberT[0].addKeyListener(this);
numberT[1].addKeyListener(this);
methodChoice.addItemListener(this);
plotChoice.addItemListener(this);
functionChoice[0].addItemListener(this);
functionChoice[1].addItemListener(this);
GridBagLayout gridbag = new GridBagLayout();//set up the layout
GridBagConstraints c = new GridBagConstraints();
setLayout(gridbag);
c.insets = new Insets(2,2,2,2); //top, left, bottom, right
int y = 1;
c.anchor = GridBagConstraints.WEST;
add(updateB, c);
c.gridx = 1;
add(imageLabel, c);
c.gridwidth = 2;
c.gridy = y++; c.gridx = 0;
add(statusLabel, c);
c.gridy = y++; c.gridx = 0;
add(methodChoice, c);
c.gridwidth = 1;
c.gridy = y++; c.gridx = 0; c.anchor = GridBagConstraints.EAST;
add(new Label("Periodicity:"), c);
c.gridx = 1; c.anchor = GridBagConstraints.WEST;
add(periodLabel, c);
for (int j=0; j<functionChoice.length; j++) {
c.gridy = y++; c.gridx = 0; c.anchor =
GridBagConstraints.EAST;
add(functionChoice[j], c);
c.gridx=1; c.anchor = GridBagConstraints.WEST;
add(numberT[j], c);
}
c.anchor = GridBagConstraints.CENTER;
c.gridy = y++; c.gridx = 0; c.anchor =
GridBagConstraints.EAST;
add(new Label("Result:"), c);
c.gridx=1; c.anchor = GridBagConstraints.WEST;
add(resultLabel, c);
c.gridwidth = 2;
c.gridy = y++; c.gridx = 0;c.anchor = GridBagConstraints.WEST;
add(toResultB, c);
c.insets = new Insets(15,2,8,2);
c.gridy = y++; c.gridx = 0; c.anchor = GridBagConstraints.EAST;
add(plotChoice, c);
this.pack();
GUI.center(this);
setResizable(false);
show();
toFront();
newImageListeners(); // we start to
work now
}
public void close() {
instance = null; //from now on,
the constructor will open a new panel
removeImageListeners();
closePlotImage(); //also
terminates the background thread
yProfile = null; //free space
yCorr = null;
super.close();
}
// these listeners are activated if the selection is changed in
the corresponding ImagePlus
public void mousePressed(MouseEvent e) { updateProfile(); }
public void mouseDragged(MouseEvent e) { updateProfile(); }
public void mouseClicked(MouseEvent e) { updateProfile(); }
// unused listeners concering actions in the corresponding
ImagePlus
public void mouseReleased(MouseEvent e) {}
public void mouseExited(MouseEvent e) {}
public void mouseEntered(MouseEvent e) {}
public void mouseMoved(MouseEvent e) {}
public void keyPressed(KeyEvent e) {}
public void keyTyped(KeyEvent e) {}
// listener for both the input fields of this panel and for the
corresponding ImagePlus
public void keyReleased(KeyEvent e) {
if (e.getSource() instanceof ImageCanvas) //nudged the
selection?
updateProfile();
else if (e.getSource() instanceof TextField)//changed a number?
updateMath();
}
// listener of the Buttons
public void actionPerformed(ActionEvent e) {
if (e.getSource() instanceof Button) {
Button b = (Button)e.getSource();
if (b == updateB) {
if (imp !=WindowManager.getCurrentImage())
newImageListeners();
else
updateProfile();
} else if (b == toResultB) {
addToResults();
}
} else {
updateMath();
}
}
//listener for the Choice items
public void itemStateChanged(ItemEvent e) {
Choice c = (Choice)e.getSource();
if (c == plotChoice) {
showPlot = c.getSelectedIndex(); // Type of Plot
changed
synchronized (this) {
if (showPlot == 0 && plotImage!=null) closePlotImage();
}
updateProfile(); // this will
open a plot image if requested
} else if (c == methodChoice) {
updateProfile();
} else updateMath(); // All others
concern the math
}
public void run() { // the
background thread for plotting. Notified to terminate by setting its
name to "...Done"
while (instance!=null && Thread.currentThread().getName
().indexOf("Done")<0) {
updatePlot();
synchronized(this) {
try {wait();}
catch(InterruptedException e) {IJ.write
("INTERRUPTED");}
}
}
}
private synchronized void makePlotImage() { //create a new
plot and start the background thread
if (plotImage==null && imp!=null && yProfile!=null) {
plotImage = new ImagePlus("Measure Periodicity Plot",
makePlotProcessor());
plotImage.show();
IJ.wait(50);
positionPlotWindow();
// start backgound thread for plotting
if (bgThread != null) { IJ.error("Measure Periodicity
internal error - duplicate thread"); return; }
bgThread = new Thread(this, "Measure Periodicity Plot");
bgThread.start();
bgThread.setPriority(bgThread.getPriority()-4);
}
}
private synchronized void closePlotImage() { //close the plot
window and terminate the background thread
if (plotImage!=null)
plotImage.getWindow().close();
if (bgThread != null) {
bgThread.setName("Measure Periodicity Plot - Done");
notifyAll(); //"Done" in the
name makes the thread terminate
bgThread = null;
plotImage = null;
}
}
private synchronized void newImageListeners() { // add listeners
for foreground image (remove old)
removeImageListeners();
imp = WindowManager.getCurrentImage(); // get the
current foreground window
if (imp == null) {
statusLabel.setText("no image");
periodicity = 0;
updateMath();
return;
}
ImageWindow win = imp.getWindow();
ImageCanvas canvas = win.getCanvas();
canvas.addMouseListener(this);
canvas.addMouseMotionListener(this);
canvas.addKeyListener(this);
lastImp = imp;
imageLabel.setText(imp.getTitle());
if (plotImage != null) positionPlotWindow();
updateProfile(); // if we are
successful, start working
}
private void removeImageListeners() {
if (lastImp == null) return;
ImageWindow win = lastImp.getWindow();
ImageCanvas canvas = win.getCanvas();
canvas.removeMouseListener(this);
canvas.removeMouseMotionListener(this);
canvas.removeKeyListener(this);
imageLabel.setText("No Image");
lastImp = null;
}
void positionPlotWindow() {
IJ.wait(500);
if (plotImage==null || imp==null) return;
ImageWindow pwin = plotImage.getWindow();
ImageWindow iwin = imp.getWindow();
if (pwin==null || iwin==null) return;
Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
Dimension plotSize = pwin.getSize();
Dimension imageSize = iwin.getSize();
if (plotSize.width==0 || imageSize.width==0) return;
Point imageLoc = iwin.getLocation();
int x = imageLoc.x+imageSize.width+10;
if (x+plotSize.width>screen.width)
x = screen.width-plotSize.width;
pwin.setLocation(x, imageLoc.y);
iwin.toFront();
}
synchronized void updateProfile() { //get line
profile and analyze it
if (!isSelection()) {
statusLabel.setText("No Line Selection");
} else {
statusLabel.setText("");
ProfilePlot profile = new ProfilePlot(imp);
yProfile = profile.getProfile();
if (yProfile==null || yProfile.length<5)
return;
xUnit = "pixels"; //the following
code is mainly for x calibration
xInc = 1;
Calibration cal = imp.getCalibration();
Roi roi = imp.getRoi();
Line line = (Line)roi;
if (cal == null) {
angle = -Math.atan2(line.y2-line.y1, line.x2-line.x1)
*180/Math.PI;
} else {
double dx = cal.pixelWidth*(line.x2 - line.x1);
double dy = cal.pixelHeight*(line.y2 - line.y1);
double length = Math.sqrt(dx*dx + dy*dy);
angle = -Math.atan2(dy, dx)*180/Math.PI;
if (yProfile.length>1) {
xInc = length/(yProfile.length-1);
xUnit = cal.getUnit();
}
}
angle -= 180*Math.round(angle/180); //reduce angles
to -90...+90 deg
yCorr = makeAutoCorrelation(yProfile); //get the
autocorrelation and the periodicity therefrom
periodicity = xInc*getPeriodicity(yCorr);
updateMath();
if (showPlot!=0) {
if (plotImage==null) {
makePlotImage();
} else {
if (plotImage.getWindow().isClosed()) {
showPlot = 0; //Plot Window
closed by user? then dispose
plotChoice.select(0);
closePlotImage();
} else {
notify(); //trigger plot
in background thread
}
}
}
}
}
void updatePlot() { //update an
existing plot
ImageProcessor ip = makePlotProcessor();
plotImage.setProcessor(null, ip);
}
ImageProcessor makePlotProcessor() { //create the
actual plot
double [] y;
double periodicity; // make local
copies of the global variables to get a consistent set
double xInc;
String title;
synchronized(this) {
periodicity = this.periodicity;
xInc = this.xInc;
if (showPlot == PROFILE)
y = (double[])yProfile.clone();
else
y = (double[])yCorr.clone();
}
int n = y.length;
double[] x = new double[n];
for (int i=0; i<n; i++)
x[i] = i*xInc;
String xLabel = "Distance (" + xUnit + ")";
String yLabel = showPlot==PROFILE ? "Height":"Autocorrelation";
Plot plot = new Plot("profile", xLabel, yLabel, x, y,
Plot.Y_GRID+Plot.Y_NUMBERS+Plot.X_TICKS+Plot.X_NUMBERS);
double[] yLine = Tools.getMinMax(y);
double[] xLine = new double[2];
if (periodicity>2*xInc) {
plot.setColor(Color.red);
double offset = showPlot==PROFILE ? xInc*getMaxOffset(y,
periodicity/xInc) : 0.;
for (double xL=offset; xL<(n-1)*xInc; xL+=periodicity) {
xLine[0] = xLine[1] = xL;
plot.addPoints(xLine,yLine,Plot.LINE);
}
plot.setColor(Color.black);
}
return plot.getProcessor();
}
void updateMath() { // do the
calculations based on the periodicity and update the panel
if (periodicity == 0) {
result = 0;
periodLabel.setText("");
resultLabel.setText("");
} else {
int digits = 4 - (int)(Math.log(periodicity)/Math.log(10));
if (digits<0) digits = 0;
periodLabel.setText(IJ.d2s(periodicity, digits)+" "+xUnit);
double temp = periodicity;
for (int j=0; j<functionChoice.length; j++) {
double num = Double.parseDouble(numberT[j].getText());
int fun = functionChoice[j].getSelectedIndex();
switch (fun) {
case 0: temp*=num; break;
case 1: temp/=num; break;
case 2: temp*=Math.sqrt(num); break;
case 3: temp/=Math.sqrt(num); break;
}
//IJ.write("math gets "+(float)temp);
}
result = temp;
if (Double.isNaN(result) || Double.isInfinite(result)) {
resultLabel.setText("");
} else {
digits = 4 - (int)(Math.log(result)/Math.log(10));
if (digits<0) digits = 0;
resultLabel.setText(IJ.d2s(result, digits)+" "+xUnit);
}
}
}
void addToResults() {
if (result == 0) return;
ResultsTable rt = ResultsTable.getResultsTable();
rt.incrementCounter();
rt.setValue("Periodicity", rt.getCounter()-1, result);
rt.setValue("Angle", rt.getCounter()-1, angle);
rt.show("Results");
}
// returns true if there is a simple line selection
boolean isSelection() {
if (imp==null)
return false;
Roi roi = imp.getRoi();
if (roi==null)
return false;
return roi.getType()==Roi.LINE;
}
double[] makeAutoCorrelation(double[] y) {
int outLength = (int)(0.8 * y.length);
double[] correlation = new double[outLength];
double avg = 0;
for (int i=0; i<y.length; i++) avg += y[i];
avg /= y.length;
double norm = 1;
for (int i=0; i<outLength; i++) {
double sum = 0;
for (int j=0; j<y.length-i; j++)
sum += (y[j]-avg)*(y[j+i]-avg);
if (i==0) norm = sum/y.length;
correlation[i] = sum/norm/(y.length-i);
}
return correlation;
}
/** Get the periodicity from an autocorrelation
* @param y the autocorrelation function
* @return periodicity (distance between high maxima in y),
0 if none found
*/
private double getPeriodicity(double[] y) {
boolean shortest = methodChoice.getSelectedItem().indexOf
("Shortest")>=0;
double period;
int bestI = 0; //index of the
best maximum
int firstI = 0; //index of the
first maximum
double bestY = 0;
double dampenHigh = 1; //factor
suppressing maxima with high periodicity
for (int i=1; i<y.length-1; i++) { //find maxima of
the autocorrelation
dampenHigh = shortest ? dampenHigh*0.9 : (1 - i*0.2/
y.length);
if (y[i]>y[i-1] && y[i]>=y[i+1]){
if (firstI==0) firstI = i;
if (y[i]*dampenHigh > bestY) { //the best
maximum so far?
bestY = y[i];
bestI = i;
}
}
}
if (bestI == 0 || bestI > 0.8*y.length) return 0;
period = bestI;
double sumPeriod = 0;
int range = firstI<10 ? 1:2; //for closely
spaced maxima, fit only a short range around the maximum
for (int i=1; period*i<y.length-2; i++) {
double x = getNearestMax(y, i*period, range);
if (Double.isNaN(x)) return 0;
sumPeriod += x; //slope is x/i;
weight proportional to i
period = sumPeriod*2/i/(i+1); //a new (and
better) estimate. Weight = sum(i) = i*(i+1)/2
if (period<2.1) return 0;
}
return period;
}
double getNearestMax(double[] y, double startX, int range) {//
find the exact point of the maximum of a function y
LineFit myFit = new LineFit(); //linear fit of the
slope, ther maximum is at zero slope
int startI = (int)Math.round(startX);
for (int i=Math.max(startI-range,0); i<Math.min(startI
+range,y.length-1); i++)
myFit.addPoint(i-startI+0.5, y[i+1]-y[i]);
double slope=myFit.getSlope();
if (slope>=0) return Double.NaN; // no minimum found
double dx = -myFit.getOffset()/slope;
if (Math.abs(dx)>2) return Double.NaN; // start point
too far from maximum
return startI + dx;
}
int getMaxOffset(double[] y, double period) { //get the
index of the maximum of y when folded (summed periodically)
int maxOffset = 0;
double maxSum = -Double.MAX_VALUE;
int repetitions = (int)Math.floor((y.length-1)/period);
//IJ.write("getMaxOffset: period="+(float)period+"
len="+y.length+" imax="+Math.floor(period)+" jmax="+repetitions);
if(repetitions<0||repetitions>y.length) {IJ.beep();return
0;} //actually an error that should not occur
for (int i=0; i<=Math.floor(period); i++) {
double sum=0;
for (int j=0; j<repetitions; j++)
sum += y[i + (int)Math.round(j*period)];
if (sum>maxSum) {
maxOffset = i;
maxSum = sum;
}
}
return maxOffset;
}
}