# 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?

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");
``````

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,

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:

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.

Hereâ€™s a â€śshortâ€ť 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("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]");
setThreshold(2, 2);
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);

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;
}
}
}
}

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
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?