Skip to content

Commit 2b5cf40

Browse files
committed
Extend path profile to timelapses
see https://forum.image.sc/t/intensity-profiles-of-traced-dendrites-over-time/121489 While at it fix use calibrated units choice not being respected all the time
1 parent cf18fe7 commit 2b5cf40

8 files changed

Lines changed: 1230 additions & 598 deletions

File tree

src/main/java/sc/fiji/snt/PathManagerUI.java

Lines changed: 91 additions & 51 deletions
Large diffs are not rendered by default.

src/main/java/sc/fiji/snt/analysis/PathProfiler.java

Lines changed: 215 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -23,15 +23,11 @@
2323
package sc.fiji.snt.analysis;
2424

2525
import java.awt.Color;
26-
import java.util.ArrayList;
27-
import java.util.Arrays;
28-
import java.util.Collections;
29-
import java.util.HashMap;
30-
import java.util.List;
31-
import java.util.Map;
26+
import java.util.*;
3227
import java.util.stream.Collectors;
3328
import java.util.stream.IntStream;
3429

30+
import net.imagej.axis.Axes;
3531
import org.apache.commons.math3.stat.descriptive.SummaryStatistics;
3632
import org.scijava.ItemVisibility;
3733
import org.scijava.command.Command;
@@ -64,6 +60,7 @@
6460
import sc.fiji.snt.util.ImgUtils;
6561
import sc.fiji.snt.util.PointInImage;
6662
import sc.fiji.snt.util.SNTColor;
63+
import sc.fiji.snt.util.TreeUtils;
6764

6865
/**
6966
* Command to retrieve Path profiles (plots of voxel intensities values along a
@@ -134,13 +131,13 @@ public class PathProfiler extends CommonDynamicCmd {
134131
private String channelString;
135132

136133
@Parameter(required = false, label = "Spatially calibrated distances")
137-
private boolean usePhysicalUnits;
134+
protected boolean usePhysicalUnits;
138135

139136
@Parameter(label = "tree")
140-
private Tree tree;
137+
protected Tree tree;
141138

142139
@Parameter(label = "dataset")
143-
private Dataset dataset;
140+
protected Dataset dataset;
144141

145142
@Parameter(label = "imp", required = false)
146143
private ImagePlus imp;
@@ -171,7 +168,7 @@ public PathProfiler(final Tree tree, final ImagePlus imp) {
171168
this.tree = tree;
172169
initContextAsNeeded();
173170
dataset = getDatasetFromImp(imp);
174-
setUnitAndAvgSept();
171+
setSpatialUnitAndAvgSept();
175172
setRadius(0);
176173
setShape(ProfileProcessor.Shape.LINE);
177174
setMetric(ProfileProcessor.Metric.MEAN);
@@ -201,7 +198,7 @@ public PathProfiler(final Tree tree, final Dataset dataset) {
201198
throw new IllegalArgumentException("Tree cannot be null");
202199
this.tree = tree;
203200
this.dataset = dataset;
204-
setUnitAndAvgSept();
201+
setSpatialUnitAndAvgSept();
205202
setRadius(0);
206203
setShape(ProfileProcessor.Shape.LINE);
207204
setMetric(ProfileProcessor.Metric.MEAN);
@@ -223,20 +220,20 @@ private void updateMsg() {
223220
if (avgSep == 0) {
224221
msg = "";
225222
} else {
226-
if (radius == 0)
223+
if (radius <= 0)
227224
msg = "Radius set by path radii";
228225
else
229-
msg = "Aprox.: " + SNTUtils.formatDouble(avgSep * radius, 2) + unit;
226+
msg = String.format("Aprox.: %.2f%s", (avgSep * radius), unit);
230227
}
231228
}
232229

233230
@SuppressWarnings("unused")
234-
private void init() {
231+
protected void init() {
235232
initContextAsNeeded();
236233
super.init(false);
237234
if (dataset == null)
238235
dataset = getDatasetFromImp(imp);
239-
setUnitAndAvgSept();
236+
setSpatialUnitAndAvgSept();
240237
try {
241238
// adjust shapeStr options
242239
final MutableModuleItem<String> mi = getInfo().getMutableInput("shapeStr", String.class);
@@ -269,12 +266,16 @@ private Dataset getDatasetFromImp(final ImagePlus imp) {
269266
return dataset;
270267
}
271268

272-
private void setUnitAndAvgSept() {
273-
unit = dataset.axis(0).unit();
269+
private void setSpatialUnitAndAvgSept() {
270+
unit = dataset.axis(dataset.dimensionIndex(Axes.X)).unit();
274271
avgSep = 0;
275-
for (int i = 0; i < dataset.numDimensions(); i++)
272+
int nSpatialAxes = 0;
273+
for (int i = 0; i < dataset.numDimensions(); i++) {
274+
if (!dataset.axis(i).type().isSpatial()) continue;
276275
avgSep += dataset.axis(i).calibratedValue(1);
277-
avgSep /= dataset.numDimensions();
276+
nSpatialAxes++;
277+
}
278+
avgSep /= nSpatialAxes;
278279
}
279280

280281
private void initContextAsNeeded() {
@@ -298,7 +299,7 @@ private void metricStrChanged() {
298299
}
299300
}
300301

301-
private void evalParameters() {
302+
protected void evalParameters() {
302303
switch (metricStr.toLowerCase()) {
303304
case "sum" -> metric = ProfileProcessor.Metric.SUM;
304305
case "min" -> metric = ProfileProcessor.Metric.MIN;
@@ -387,7 +388,7 @@ public void run() {
387388
}
388389
}
389390

390-
private List<Integer> getChannels() {
391+
protected List<Integer> getChannels() {
391392
if (channelString == null || channelString.trim().isEmpty() || "all".equalsIgnoreCase(channelString.trim()))
392393
return getAllChannels();
393394
final List<String> stringChannels = new ArrayList<>(Arrays.asList(channelString.split("\\s*([,\\s])\\s*")));
@@ -484,8 +485,24 @@ public void assignValues(final Path p) throws IllegalArgumentException {
484485
* @throws IllegalArgumentException if image does not contain the path's channel
485486
*/
486487
public <T extends RealType<T>> void assignValues(final Path p, final int channel) throws ArrayIndexOutOfBoundsException {
488+
assignValues(p, channel, p.getFrame()-1);
489+
//ImgUtils.getCtSlice3d(dataset, channel, channel);
490+
}
491+
492+
/**
493+
* Retrieves pixel intensities at each node of the Path storing them as Path
494+
* {@code values}
495+
*
496+
* @param channel the channel to be parsed (base-0 index)
497+
* @param frame the frame to be parsed (base-0 index)
498+
* @param p the Path to be profiled
499+
* @see Path#setNodeValues(double[])
500+
*
501+
* @throws IllegalArgumentException if image does not contain the path's channel
502+
*/
503+
public <T extends RealType<T>> void assignValues(final Path p, final int channel, final int frame) throws ArrayIndexOutOfBoundsException {
487504
validateChannelRange(channel);
488-
final RandomAccessibleInterval<T> rai = ImgUtils.getCtSlice(dataset, channel, p.getFrame() - 1);
505+
final RandomAccessibleInterval<T> rai = ImgUtils.getCtSlice(dataset, channel, frame);
489506
final ProfileProcessor<T> processor = new ProfileProcessor<>(rai, p);
490507
processor.setShape(shape);
491508
processor.setRadius(radius);
@@ -578,14 +595,30 @@ public Map<String, List<Double>> getValues(final Path p) {
578595
/**
579596
* Gets the profile for the specified path as a map of lists, with distances (or
580597
* indices) stored under {@link #X_VALUES} ({@value #X_VALUES}) and intensities
581-
* under {@link #Y_VALUES} ({@value #Y_VALUES}).
598+
* under {@link #Y_VALUES} ({@value #Y_VALUES}). If dataset is a time-lapse, the
599+
* path's assigned frame is profiled.
582600
*
583601
* @param p the path to be profiled
584602
* @param channel the channel to be parsed (base-0 index)
585603
* @return the profile map
586604
*/
587605
public Map<String, List<Double>> getValues(final Path p, final int channel) {
588-
if (!p.hasNodeValues()) assignValues(p, channel);
606+
return getValues(p, channel, p.getFrame()-1);
607+
}
608+
609+
/**
610+
* Gets the profile for the specified path as a map of lists, with distances (or
611+
* indices) stored under {@link #X_VALUES} ({@value #X_VALUES}) and intensities
612+
* under {@link #Y_VALUES} ({@value #Y_VALUES}).
613+
*
614+
* @param p the path to be profiled
615+
* @param channel the channel to be parsed (base-0 index)
616+
* @param frame the frame to be parsed (base-0 index)
617+
* @return the profile map
618+
*/
619+
public Map<String, List<Double>> getValues(final Path p, final int channel, final int frame) {
620+
621+
if (!p.hasNodeValues()) assignValues(p, channel, frame);
589622
final List<Double> xList = new ArrayList<>();
590623
final List<Double> yList = new ArrayList<>();
591624

@@ -664,12 +697,12 @@ private ColorRGB[] getSeriesColorsRGB() {
664697
return colors;
665698
}
666699

667-
private String getXAxisLabel() {
700+
protected String getXAxisLabel() {
668701
return (nodeIndices) ? "Node indices"
669702
: String.format("Distance (%s)", tree.getProperties().getProperty(Tree.KEY_SPATIAL_UNIT, "? units"));
670703
}
671704

672-
private String getYAxisLabel(final int channel) {
705+
protected String getYAxisLabel(final int channel) {
673706
final boolean detailed = shape != Shape.NONE;
674707
final StringBuilder sb = new StringBuilder();
675708
if (channel > 0 && dataset.getChannels() > 1) {
@@ -678,7 +711,13 @@ private String getYAxisLabel(final int channel) {
678711
sb.append(dataset.getValidBits()).append("-bit ");
679712
if (detailed) {
680713
sb.append("Int. (").append(metric).append("; ");
681-
sb.append(shape).append(", r=").append((radius==0) ? "Node radius" : radius + "px");
714+
sb.append(shape).append(", r=");
715+
if (radius <=0) {
716+
sb.append("Node radius");
717+
} else {
718+
sb.append(String.format("%.2f", (usePhysicalUnits) ? (avgSep * radius) : radius))
719+
.append((usePhysicalUnits) ? unit : "px");
720+
}
682721
sb.append(")");
683722
} else {
684723
sb.append("Intensity");
@@ -826,6 +865,155 @@ private boolean treeIsColorMapped() {
826865
return false;
827866
}
828867

868+
/**
869+
* Resamples the intensity profile of a path onto a uniform distance grid via
870+
* linear interpolation. Entries beyond the path's actual length are set to
871+
* {@link Double#NaN}, which allows variable-length paths to be compared in
872+
* the same matrix without padding with zeros.
873+
*
874+
* @param p the path to profile (values are assigned if not yet set)
875+
* @param channel the channel to sample (base-0 index)
876+
* @param nSamples the number of output grid points (&gt;= 2)
877+
* @param gridMax the upper bound of the distance grid. Use the path's actual
878+
* length for a per-path-normalized grid (0..length), or a
879+
* shared maximum across all paths for an absolute grid
880+
* @return array of length {@code nSamples} with interpolated intensities
881+
*/
882+
public double[] getResampledValues(final Path p, final int channel,
883+
final int nSamples, final double gridMax) {
884+
if (nSamples < 2)
885+
throw new IllegalArgumentException("nSamples must be >= 2");
886+
assignValues(p, channel); // uses p.getFrame()-1
887+
return resampleAssignedValues(getValues(p, channel), nSamples, gridMax);
888+
}
889+
890+
/**
891+
* Resamples the intensity profile of a path onto a uniform distance grid via
892+
* linear interpolation. Entries beyond the path's actual length are set to
893+
* {@link Double#NaN}, which allows variable-length paths to be compared in
894+
* the same matrix without padding with zeros.
895+
*
896+
* @param p the path to profile (values are assigned if not yet set)
897+
* @param channel the channel to sample (base-0 index)
898+
* @param frame the frame to sample (base-0 index)
899+
* @param nSamples the number of output grid points (&gt;= 2)
900+
* @param gridMax the upper bound of the distance grid. Use the path's actual
901+
* length for a per-path-normalized grid (0..length), or a
902+
* shared maximum across all paths for an absolute grid
903+
* @return array of length {@code nSamples} with interpolated intensities
904+
*/
905+
public double[] getResampledValues(final Path p, final int channel, final int frame,
906+
final int nSamples, final double gridMax) {
907+
if (nSamples < 2)
908+
throw new IllegalArgumentException("nSamples must be >= 2");
909+
assignValues(p, channel, frame);
910+
return resampleAssignedValues(getValues(p, channel, frame), nSamples, gridMax);
911+
}
912+
913+
private double[] resampleAssignedValues(final Map<String, List<Double>> values,
914+
final int nSamples, final double gridMax) {
915+
if (nSamples < 2)
916+
throw new IllegalArgumentException("nSamples must be >= 2");
917+
// Always assign values explicitly: PointInImage.v defaults to 0.0, not NaN,
918+
// so hasNodeValues() returns true for unsampled paths, causing getValues()
919+
// to silently skip assignValues() and return all-zero intensities.
920+
final List<Double> srcX = values.get(X_VALUES);
921+
final List<Double> srcY = values.get(Y_VALUES);
922+
final double pathLen = srcX.getLast();
923+
final double step = gridMax / (nSamples - 1);
924+
final double[] result = new double[nSamples];
925+
int j = 0; // pointer into srcX/srcY
926+
for (int i = 0; i < nSamples; i++) {
927+
final double x = i * step;
928+
if (x > pathLen + 1e-9) {
929+
result[i] = Double.NaN;
930+
continue;
931+
}
932+
// advance j so srcX[j] <= x < srcX[j+1]
933+
while (j < srcX.size() - 2 && srcX.get(j + 1) <= x) j++;
934+
if (j >= srcX.size() - 1) {
935+
result[i] = srcY.get(srcX.size() - 1);
936+
} else {
937+
final double x0 = srcX.get(j), x1 = srcX.get(j + 1);
938+
final double y0 = srcY.get(j), y1 = srcY.get(j + 1);
939+
result[i] = (x1 == x0) ? y0 : y0 + (y1 - y0) * (x - x0) / (x1 - x0);
940+
}
941+
}
942+
return result;
943+
}
944+
945+
/**
946+
* Builds a multi-frame intensity profile matrix from a list of matched paths
947+
* (one per time frame, ordered by frame number), suitable for kymograph-style
948+
* visualization via {@link SNTChart#showHeatmap}.
949+
* <p>
950+
* The returned matrix has dimensions [nPaths][nSamples], where
951+
* {@code matrix[i][j]} is the interpolated intensity at distance sample j for
952+
* the i-th path. Entries beyond a path's actual length are {@link Double#NaN}.
953+
* </p>
954+
*
955+
* @param paths ordered list of paths, one per time frame
956+
* @param channel the channel to sample (base-0 index)
957+
* @param nSamples number of distance samples (columns in the result)
958+
* @param normalizeDistance if true, distances are normalized to [0,1] so all
959+
* rows span the full width regardless of path length;
960+
* if false, absolute distances are used and the grid
961+
* max equals the longest path in the list
962+
* @return a 2D array [nPaths][nSamples] of intensity values
963+
*/
964+
public double[][] getMultiFrameProfile(final List<Path> paths, final int channel,
965+
final int nSamples, final boolean normalizeDistance) {
966+
if (paths == null || paths.isEmpty())
967+
throw new IllegalArgumentException("Paths list cannot be null or empty");
968+
// For absolute distances, share a single grid max so all rows are comparable.
969+
// For normalized distances, each path is resampled over its OWN length so
970+
// the grid always spans [0, pathLen] and maps to [0, 1] conceptually --
971+
// passing gridMax = 1.0 is wrong because source x-values are in physical
972+
// units (e.g. um), which would confine sampling to the first sub-unit segment.
973+
final double sharedGridMax;
974+
if (normalizeDistance) {
975+
sharedGridMax = -1; // unused; each path uses its own length below
976+
} else {
977+
double maxLen = 0;
978+
for (final Path p : paths) {
979+
final double len = p.getLength();
980+
if (len > maxLen) maxLen = len;
981+
}
982+
sharedGridMax = maxLen;
983+
}
984+
final double[][] matrix = new double[paths.size()][nSamples];
985+
for (int i = 0; i < paths.size(); i++) {
986+
final Path p = paths.get(i);
987+
final double gridMax = normalizeDistance ? p.getLength() : sharedGridMax;
988+
matrix[i] = getResampledValues(p, channel, nSamples, gridMax);
989+
}
990+
return matrix;
991+
}
992+
993+
/**
994+
* Gets the time profile for the specified path as a list of {@link #getValues(Path, int, int)} profiles,
995+
* with one entry per frame of the dataset.
996+
*
997+
* @param path the path to be profiled
998+
* @param channel the channel to be parsed (base-0 index)
999+
* @return the list of profile values (one entry per frame)
1000+
*/
1001+
public List<Map<String, List<Double>>> getTimeProfile(final Path path, final int channel) {
1002+
if (path == null || path.size()==0)
1003+
throw new IllegalArgumentException("Path cannot be null or empty");
1004+
final Map<Path, double[]> nodeValuesSnapshot = TreeUtils.snapshotNodeValues(new Tree(List.of(path)));
1005+
final List<Map<String, List<Double>>> result = new ArrayList<>();
1006+
try {
1007+
for (int frame = 0; frame < dataset.getFrames(); frame++) {
1008+
assignValues(path, channel, frame);
1009+
result.add(getValues(path, channel, frame));
1010+
}
1011+
} finally {
1012+
TreeUtils.restoreNodeValues(nodeValuesSnapshot);
1013+
}
1014+
return result;
1015+
}
1016+
8291017
/* IDE debug method **/
8301018
public static void main(final String[] args) {
8311019
final ImageJ ij = new ImageJ();

src/main/java/sc/fiji/snt/analysis/PathStatistics.java

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ public int getNBranches() {
143143
@Override
144144
public Number getMetric(final String metric) throws UnknownMetricException {
145145
if ("Path ID".equalsIgnoreCase(metric))
146-
return (tree.size() == 1) ? tree.list().get(0).getID() : Double.NaN;
146+
return (tree.size() == 1) ? tree.list().getFirst().getID() : Double.NaN;
147147
return super.getMetric(metric);
148148
}
149149

@@ -373,8 +373,8 @@ public void measureIndividualPaths(final Collection<String> metrics, final boole
373373
table.set(getCol(metric), row, new PathStatistics(path).getMetric(metric, path));
374374
});
375375
});
376-
if (summarize && table instanceof SNTTable) {
377-
((SNTTable) table).summarize();
376+
if (summarize && table != null) {
377+
table.summarize();
378378
}
379379
updateAndDisplayTable();
380380
}

0 commit comments

Comments
 (0)