Get XY coordinates from FindPeaks result

Hello,

is there already some implementation to get the XY coordinates in the original image, of the maximum peaks of a plot profile?

Thank you in advance.
Mafalda

Hi Mafalda!

There’s a macro function (Array.findMaxima) to get the maxima from an array, within a tolerance (there’s also one for the minima). You can combine it with getProfile to extract the profile maxima, then use the ends of the profile line to map those coordinates back to the image, like in this example:

tolerance=100;

if (selectionType()!=5) exit("A line ROI is needed");

getSelectionCoordinates(xpoints, ypoints);

profile=getProfile();
peaks=Array.findMaxima(profile, tolerance); //get the maximima within tolerance

xpeaks=newArray(peaks.length);
ypeaks=newArray(peaks.length);

for (i=0; i<peaks.length; i++){
	q=peaks[i]/(profile.length-1);
	xpeaks[i]=xpoints[0]+(xpoints[1]-xpoints[0])*(q);
	ypeaks[i]=ypoints[0]+(ypoints[1]-ypoints[0])*(q);
	}

if(peaks.length>0){
	Array.show(xpeaks, ypeaks); // show the coordinates of the peaks
	makeSelection("points", xpeaks, ypeaks); // draw a point roi with those coords
	}
else exit("No peaks found"); 

maxima

Is this the sort of solution you were looking for?

Cheers!
Nico

3 Likes

Hi Mafalda,
concerning your first question, a list of maxima of an image (not a plot), you can use “Find Maxima” and select “List” as “Output type”. Then you get a Results table with the x&y coordinates of the maxima. Use the Command Recorder (Plugins>Macros>Record) to see what the macro command looks like.
You can retrieve the table values with the macro command Table.get(columnName, rowIndex).
– Michael

Hello,

thank you all for your replies. I already manage to get what I asked for, but I’m having some unexpected results.
Here is the code:

///roi is selected
run("Interpolate", "interval=1");
run("Area to Line");
Roi.getCoordinates(xpoints, ypoints);
run("Plot Profile");
run("Find Peaks", "min._peak_amplitude=30000 min._peak_distance=50 min._value=[] max._value=[] list");
Table.rename("Plot Values.csv", "Results");

for(res =0; res<nResults; res++){
	xx =  getResult("X1", res);
	if (isNaN(xx)){
		number_fibers = res+1;
		break;
	}
}
print("Number of fibers = ", number_fibers);
positions_X = newArray(number_fibers);
positions_Y = newArray(number_fibers);
for (i = 0; i < number_fibers; i++) {
	print(i);
	positions_X[i] = xpoints[getResult("X1", i)];
	positions_Y[i] = ypoints[getResult("X1", i)]; 
	
}

I’m asking for points that are separated more than min._peak_distance=50, but I’m having closer points. Do you have any idea why?
C1-rgb-1.tif (101.9 KB)

Thank you

Hi Mafalda,

I’m a little confused regarding the goal of your macro. Could you give us some example of the type of analysis you are trying to achieve?

The second line of your macro (run("Area to Line")), gets an area Roi and makes a snip in its perimeter, effectively turning it into a line Roi. Is this what your are looking for? If you use this on a thresholded area in the image you provided (a fiber?), you would end up with an intensity measurement along the perimeter, not along its axis. Also, you would have two readings of the same section (one on each side). In fact, as you are not using the “exclude peaks on edges of plot” option in the BAR plugin, if you have a “high spot” in the place where the command splits the perimeter, you would get two adjacent points selected (as the perimeter wraps up in the same spot).

Cheers,
Nico

Hi Nico,

Thanks for your concern.

I’m trying to count the number of fibers that are crossing a line (that was previously obtained by thresholding and manipulating rois).

You are right in most of your assumptions: The Roi is not a line but an area (I didn’t manage to get only one). That’s why I need to convert area to Line.

That’s not really true since I’m only having two closed detections in some spots:

crop_roi_peaks
crop_roi_spots

Cheers
Mafalda

Hi Mafalda,

In fact, that’s actually the case. If you look closely, the profile you showed is almost a mirrored image. You are counting the crossings twice, gowing from one end of the ROI to the other, and once again on the way back.

You can also check the following: the total length of the profile (about 60000 units) should be twice its expected length.

The roi line shown in the images is in fact double, a few pixels apart. Do you use the Make band command at some point? I can’t find a simple built in solution to get this kind of thin areas into a single line. I know @superresolusian was working on something similiar a while back . Also, this post form @iarganda might provide a useful solution, by first skeletonizing the area and then getting the longest line.

I hope this helps.

Cheers!
Nico

2 Likes

Ok, I couldn’t resist myself from trying my hand at it. :grin:

Here’s a “short” :roll_eyes: macro to transform a thin area ROI into a single line:

//	Transform an area ROI to a central (skeleton) line

setBatchMode(true);

// get all the loose ends
run("Create Mask");
selectWindow("Mask");
run("Set Scale...", "distance=0 known=0"); // in case the original image has a scale set
setOption("BlackBackground", false);
run("Skeletonize");
run("Divide...", "value=255");
run("Duplicate...", "title=ends");
run("Convolve...", "text1=[1 1 1\n1 1 1\n1 1 1\n]");
imageCalculator("Multiply", "ends","Mask");
setThreshold(2, 2);
run("Convert to Mask");
run("Points from Mask");
getSelectionCoordinates(xEnds, yEnds);
close("ends");

if(xEnds.length>2) exit("This ROI produces a branched skeleton");

// trace the skeleton
getStatistics(area, mean);
np=round(area*mean);

run("Restore Selection");

xc=newArray(np);
yc=newArray(np);

xc[0]=xEnds[0];
yc[0]=yEnds[0];

dx=newArray(-1, 0, 1,-1, 1,-1, 0, 1);
dy=newArray(-1,-1,-1, 0, 0, 1, 1, 1);
selectWindow("Mask");

for (i=0; i<np; i++){
	v=newArray(8);
	last=maxOf(0, i-1);
	lastx=xc[last];
	lasty=yc[last];
	for (j=0; j<8; j++){
		x=xc[i]+dx[j];
		y=yc[i]+dy[j];
		v[j]=getPixel(x, y);
		if(v[j]>0){
			if(!(x==lastx && y==lasty)){		
				xc[i+1]=x;
				yc[i+1]=y;
				}
			}
		} 
	}

close("Mask");
setBatchMode(false);
makeSelection("freeline", xc, yc); 

This whole block should replace the run("Area to Line") command in your macro.

Let me know if it works as expected.

Cheers!
Nico

Hi Nico,

very good you couldn’t resist from trying :slight_smile:
However I’m having a branched skeleton and stoping at:

if(xEnds.length>2) exit("This ROI produces a branched skeleton");

Any suggestion?

Mafalda

That’s odd. What steps do you follow to get the thin area ROI? I feel there might be something that can be modified to avoid overcomplicating the downstream processing.

Also, can you spot where the skeleton is branching?