Path.java

/*-
 * #%L
 * Fiji distribution of ImageJ for the life sciences.
 * %%
 * Copyright (C) 2010 - 2026 Fiji developers.
 * %%
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public
 * License along with this program.  If not, see
 * <http://www.gnu.org/licenses/gpl-3.0.html>.
 * #L%
 */

package sc.fiji.snt;

import ij.ImagePlus;
import ij.measure.Calibration;
import ij3d.Content;
import ij3d.Image3DUniverse;
import ij3d.Pipe;
import ij3d.Utils;
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.stat.regression.SimpleRegression;
import org.apache.commons.math3.util.MathUtils;
import org.scijava.util.ColorRGB;
import org.scijava.util.ColorRGBA;
import org.jogamp.vecmath.Color3f;
import org.jogamp.vecmath.Point3f;
import org.jogamp.vecmath.Vector3d;
import sc.fiji.snt.analysis.PathProfiler;
import sc.fiji.snt.annotation.BrainAnnotation;
import sc.fiji.snt.gui.GuiUtils;
import sc.fiji.snt.hyperpanes.MultiDThreePanes;
import sc.fiji.snt.util.*;

import java.awt.*;
import java.util.List;
import java.util.function.DoublePredicate;
import java.util.*;

/**
 * This class represents a traced segment (i.e., a <i>Path</i>) in a
 * reconstruction. It has methods to manipulate its points (nodes) with
 * sup-pixel accuracy, including drawing them onto threePane-style canvases,
 * etc.
 * <p>
 * Paths are non-branching sequences of adjacent points (including diagonals) in
 * the image. Branches and joins are supported by attributes of paths that
 * specify that they begin on (or end on) other paths.
 * </p>
 * <p>
 * In SNT, a Path can exist in two versions or flavors: itself and a fitted
 * version generated by {@link PathFitter}. Because fitting may fail around
 * certain nodes, the fitted version may have a different number of nodes
 * relatively to its non-fitted version
 * </p>
 **/
public class Path implements Comparable<Path>, Cloneable {

	static { net.imagej.patcher.LegacyInjector.preinit(); } // required for _every_ class that imports ij. classes

	// https://swc-specification.readthedocs.io/en/latest/
	/** SWC type flag specifying {@value #SWC_UNDEFINED_LABEL} */
	public static final int SWC_UNDEFINED = 0;
	/** SWC type flag specifying {@value #SWC_SOMA_LABEL} */
	public static final int SWC_SOMA = 1;
	/** SWC type flag specifying {@value #SWC_AXON_LABEL} */
	public static final int SWC_AXON = 2;
	/** SWC type flag specifying {@value #SWC_DENDRITE_LABEL} */
	public static final int SWC_DENDRITE = 3;
	/** SWC type flag specifying {@value #SWC_APICAL_DENDRITE_LABEL} */
	public static final int SWC_APICAL_DENDRITE = 4;
	/** SWC type flag specifying {@value #SWC_CUSTOM_LABEL} */
	public static final int SWC_CUSTOM = 5;
	/** SWC type flag specifying {@value #SWC_UNSPECIFIED_LABEL} */
	public static final int SWC_UNSPECIFIED = 6;
	/** SWC type flag specifying {@value #SWC_GLIA_PROCESS_LABEL} */
	public static final int SWC_GLIA_PROCESS = 7;
	/** SWC type flag specifying {@value #SWC_CUSTOM2_LABEL} */
	public static final int SWC_CUSTOM2 = 8;

	/** Deprecated. No longer part of the SWC specification */
	@Deprecated
	public static final int SWC_FORK_POINT = -5;
	/** Deprecated. No longer part of the SWC specification */
	@Deprecated
	public static final int SWC_END_POINT = -6;
	
	/** String representation of {@link Path#SWC_UNDEFINED} */
	public static final String SWC_UNDEFINED_LABEL = "undefined";
	/** String representation of {@link Path#SWC_SOMA} */
	public static final String SWC_SOMA_LABEL = "soma";
	/** String representation of {@link Path#SWC_AXON} */
	public static final String SWC_AXON_LABEL = "axon";
	/** String representation of {@link Path#SWC_DENDRITE} */
	public static final String SWC_DENDRITE_LABEL = "(basal) dendrite";
	/** String representation of {@link Path#SWC_APICAL_DENDRITE} */
	public static final String SWC_APICAL_DENDRITE_LABEL = "apical dendrite";
	/** String representation of {@link Path#SWC_CUSTOM} */
	public static final String SWC_CUSTOM_LABEL = "custom";
	/** String representation of {@link Path#SWC_UNSPECIFIED} */
	public static final String SWC_UNSPECIFIED_LABEL = "unspecified neurite";
	/** String representation of {@link Path#SWC_GLIA_PROCESS} */
	public static final String SWC_GLIA_PROCESS_LABEL = "glia process";
	/** String representation of {@link Path#SWC_CUSTOM2} */
	public static final String SWC_CUSTOM2_LABEL = "custom (2)";

	// Default sampling interval for 3D visualization - can be overridden based on path characteristics
	protected static final int DEFAULT_SAMPLING_INTERVAL = 2;
	
	/**
	 * @deprecated Use {@link #DEFAULT_SAMPLING_INTERVAL} instead
	 */
	@Deprecated
	protected static final int noMoreThanOneEvery = DEFAULT_SAMPLING_INTERVAL;
	
	/**
	 * Calculates an appropriate sampling interval for 3D visualization based on path characteristics.
	 * Uses path length and node density to determine optimal sampling.
	 * 
	 * @return the recommended sampling interval
	 */
	 private int calculateSamplingInterval() {
		if (size() <= DEFAULT_SAMPLING_INTERVAL) return 1;
		// For short paths, use minimal sampling
		if (size() <= 10) return 1;
		// For longer paths, calculate based on average inter-node distance
		final double avgDistance = getLength() / (size() - 1);
		final double minSpacing = getMinimumSeparation();
		// If nodes are very close together, increase sampling interval
		if (avgDistance < minSpacing * 2) {
			return Math.min(size() / 10, 5); // Max interval of 5
		}
		return DEFAULT_SAMPLING_INTERVAL;
	}

	/**
	 * Enhanced PointInImage that stores all node-specific properties.
	 * This eliminates the need for parallel arrays and provides better
	 * memory efficiency and data encapsulation.
	 */
	public static class PathNode extends PointInImage implements Cloneable {
		
		private Color nodeColor;
		private double[] tangent; // [x, y, z]

		public PathNode(final double x, final double y, final double z) {
			super(x, y, z);
		}
		
		public PathNode(final SNTPoint point) {
			super(point.getX(), point.getY(), point.getZ());
			if (point instanceof PointInImage pim) {
                this.v = pim.v;
				this.setAnnotation(pim.getAnnotation());
				this.setHemisphere(pim.getHemisphere());
                this.radius = pim.radius;
            }
		}
		
		// Node color management
		public Color getColor() { return nodeColor; }
		public boolean hasColor() { return nodeColor != null; }
		public void setColor(final Color color) { this.nodeColor = color; }

		// Radius management
		public double getRadius() { return radius; }
		public boolean hasRadius() { return radius > 0d; }
		public void setRadius(final double radius) { this.radius = radius; }

		// Tangent management
		public double[] getTangent() { return tangent; }
		public boolean hasTangent() { return tangent != null; }
		public void setTangent(final double[] tangent) { this.tangent = tangent; }
		public void setTangent(final double tx, final double ty, final double tz) {
			this.tangent = new double[]{tx, ty, tz};
		}

        private Vector3d asVector() {
            return new Vector3d(x, y, z);
        }

		/**
		 * Creates a deep copy of this PathNode.
		 * <p>
		 * This method creates a complete copy of the node including all properties
		 * such as color, radius, and tangent information. The cloned node will have
		 * the same spatial coordinates and metadata as the original.
		 * </p>
		 * 
		 * @return a new PathNode that is a deep copy of this node
		 */
		@Override
		public PathNode clone() {
			try {
				// Call super.clone() to get a shallow copy
				PathNode cloned = (PathNode) super.clone();
				
				// Deep clone mutable fields
				cloned.tangent = this.tangent != null ? this.tangent.clone() : null;
				
				return cloned;
			} catch (final Exception e) {
				// This should never happen since we implement Cloneable
				throw new AssertionError("Clone not supported", e);
			}
		}
	}

	/* Path properties */
	// Single collection of enhanced nodes (replaces all parallel arrays previous to SNTv5)
	private final List<PathNode> nodes;
	/*
	 * Path identifiers: this Path's id is stored in (lower) bits 15-0. Tree id in
	 * the (upper) bits 31-16. NB: should only be assigned by PathAndFillManager.
	 * NB: A regular path should have a path ID > 0, because if present, its fitted
	 * version will have an ID of -path ID. see {@link #setFitted(Path)}
	 */
	private long id = -1L;
	private String treeLabel;

	// NB: The legacy 3D viewer requires always a unique name
	private String name;
	// Path based ordering akin of reverse Horton-Strahler numbers
	private int order = 1;
	// The SWC-type flag of this path
	int swcType = SWC_UNDEFINED;
	// is this path selected in the UI?
	private boolean selected;
	// the node being edited when in 'Analysis mode'
	private int editableNodeIndex = -1;
	// whether node being edited is locked
	private boolean editableNodeLocked;

	// the display offset for this Path in a tracing canvas. This allows to overlay the path
	// away from the signal without changing its coordinates
	protected PointInCanvas canvasOffset = new PointInCanvas(0, 0, 0);

	// the channel and frame associated with this Path (1-based indices)
	private final int[] ctPosition;
	// the number of spines or varicosities associated with this path.
	//TODO: This should be move to PathNode so that such count can be assigned to individual nodes
	private int spinesOrVaricosities;

	/* Spatial calibration definitions */
	protected double x_spacing;
	protected double y_spacing;
	protected double z_spacing;
	protected String spacing_units;

	/* Branching */
	protected Path parentPath;
	protected int branchPointIndex = -1;
	// This is a symmetrical relationship, showing
	// all the other paths this one is joined to...
	protected ArrayList<Path> connectedPaths;
	// We almost always impose tree structure on the Path graph.
	// When this is done, we regenerated this list. This should
	// always be a subset of 'somehowJoins'...
	protected ArrayList<Path> children;

	/* Fitting (Path refinement) */
	// If this path has a fitted version, this is it
	protected Path fitted;
	// Prefer the fitted flavor of this path
	private boolean useFitted;
	// If this path is a fitted version of another one, this is the original
	protected Path fittedVersionOf;

	/* General color definitions (in addition to node colors) */
	private Color color;
	private Color3f realColor;
	private boolean hasCustomColor;

	private final List<PathChangeListener> changeListeners;

	/**
	 * Instantiates a new path under default settings (isotropic 1um pixel spacing)
	 */
	public Path() {
		this(1, 1, 1, GuiUtils.micrometer());
	}

	/**
	 * Instantiates a new path.
	 *
	 * @param x_spacing Pixel width in spacing_units
	 * @param y_spacing Pixel height in spacing_units
	 * @param z_spacing Pixel depth in spacing_units
	 * @param spacing_units the length unit in physical world units (typically "um").
	 */
	public Path(final double x_spacing, final double y_spacing,
		final double z_spacing, final String spacing_units) {
		this.x_spacing = x_spacing;
		this.y_spacing = y_spacing;
		this.z_spacing = z_spacing;
		this.spacing_units = SNTUtils.getSanitizedUnit(spacing_units);
		this.nodes = new ArrayList<>();
		connectedPaths = new ArrayList<>();
		children = new ArrayList<>();
		ctPosition = new int[] {1, 1};
		changeListeners = new ArrayList<>();
	}

	/* (non-Javadoc)
	 * @see java.lang.Comparable#compareTo(java.lang.Object)
	 */
	@Override
	public int compareTo(final Path o) {
		return Long.compare(id, o.id);
	}

	/**
	 * Gets the identifier of this Path
	 *
	 * @return the identifier
	 */
	public int getID() {
		return (int) id; // (lower) bits 15-0
	}

	/**
	 * Gets the identifier of the {@link Tree} associated to this Path (if any).
	 *
	 * @return the Tree identifier
	 */
	public int getTreeID() {
		return (int) (id >> 32); // (upper) bits 31-16
	}

	protected void setIDs(final int pathID, final int treeID) {
		final int oldId = getID();
		// pathID to fill lower bits, treeID to fill upper bits
		id = (((long) treeID) << 32) | (pathID & 0xffffffffL);
		changeListeners.forEach(l -> l.pathChanged(
				new PathChangeEvent(this, PathChangeEvent.EventType.ID_CHANGED, oldId, getID())));
	}

	protected void setTreeLabel(final String treeLabel) {
		this.treeLabel = treeLabel;
		if (getFitted() != null) getFitted().setTreeLabel(treeLabel); 
	}

	public String getTreeLabel() {
		if (isFittedVersionOfAnotherPath()) {
			return fittedVersionOf.getTreeLabel();
		}
		return (treeLabel == null) ? String.format("Arbor %d C%dT%d", getTreeID(), getChannel(), getFrame()) : treeLabel;
	}

	/**
	 * Specifies a translation offset when rendering this Path in a
	 * {@link TracerCanvas}. Path coordinates remain unaltered.
	 *
	 * @param canvasOffset the x,y,z coordinates (pixel-based) specifying the
	 *          translation offset
	 */
	public void setCanvasOffset(final PointInCanvas canvasOffset) {
		this.canvasOffset = canvasOffset;
		if (getFitted() != null) getFitted().setCanvasOffset(canvasOffset);
	}

	/**
	 * Returns the translation offset used to render this Path in a
	 * {@link TracerCanvas}.
	 *
	 * @return the rendering offset (in pixel coordinates)
	 */
	public PointInCanvas getCanvasOffset() {
		canvasOffset.onPath = this;
		return canvasOffset;
	}

	/**
	 * Gets the parent path that this path branches from.
	 * <p>
	 * Returns the path from which this path originates as a branch. For primary
	 * (root) paths, this method returns null since they don't branch from any
	 * other path. The parent-child relationship is established when a path is
	 * created as a branch using {@link #setBranchFrom(Path, PointInImage)}.
	 * </p>
	 * 
	 * @return the parent path, or null if this is a primary (root) path
	 * 
	 * @see #setBranchFrom(Path, PointInImage) for establishing parent-child relationships
	 * @see #isPrimary() for checking if this is a root path
	 * @see #getBranchPoint() for the specific branching point
	 */
	public Path getParentPath() {
		return parentPath;
	}

	/**
	 * @deprecated  use {@link #getParentPath()} instead
	 */
	@Deprecated
	public Path getStartJoins() {
		return getParentPath();
	}


	/**
	 * Gets the branch point (junction) where this path starts, i.e., connects to its parent path.
	 * 
	 * @return the start junction point, or null if this is a primary path
	 * @see #getBranchPoints()
	 */
	public PointInImage getBranchPoint() {
		if (parentPath == null || branchPointIndex < 0 || branchPointIndex >= parentPath.size()) {
			return null;
		}
		return parentPath.getNode(branchPointIndex);
	}
	
	/**
	 * Gets the index of the branch point in the parent path.
	 * 
	 * @return the branch point index, or -1 if this is a primary path
	 */
	public int getBranchPointIndex() {
		return branchPointIndex;
	}
	
	/**
	 * @deprecated Use {@link #getBranchPoint()} instead.
	 *
	 * <p>
	 * NB: Earlier versions of Simple Neurite Tracer allowed paths to connect at either the start
	 * point (getStartJoinsPoint()) or end-point (getStartJoinsPoint()). Useful for some
	 * applications but a huge headache to handle.
	 * </p>
	 */
	@Deprecated
	public PointInImage getStartJoinsPoint() {
		return getBranchPoint();
	}

	/**
	 * Sets this Path's name. Set it to null or {@code ""}, to reset it to the
	 * default.
	 *
	 * @param newName the new name.
	 * @see #getName()
	 */
	public void setName(final String newName) {
		final String oldName =  getName();
		this.name = newName;
		getName(); // assign default if newName is null
		changeListeners.forEach(l -> l.pathChanged(
				new PathChangeEvent(this, PathChangeEvent.EventType.NAME_CHANGED, oldName, getName())));
	}

	/**
	 * Gets this Path's name.
	 *
	 * @return the name. If no name as been set, the default name is returned.
	 * @see #setName(String)
	 */
	public String getName() {
		if (name == null || name.isEmpty()) name = "Path " + getID();
		return name;
	}

	protected static String pathsToIDListString(final ArrayList<Path> a) {
		final StringBuilder s = new StringBuilder();
		final int n = a.size();
		for (int i = 0; i < n; ++i) {
			s.append(a.get(i).getID());
			if (i < n - 1) {
				s.append(",");
			}
		}
		return s.toString();
	}

	protected String somehowJoinsAsString() {
		return pathsToIDListString(connectedPaths);
	}

	protected String childrenAsString() {
		return pathsToIDListString(children);
	}

	protected void setChildren(final Set<Path> pathsLeft) {
		// Set the children of this path in a breadth first fashion:
		children.clear();
		for (final Path c : connectedPaths) {
			if (pathsLeft.contains(c)) {
				children.add(c);
				pathsLeft.remove(c);
			}
		}
		for (final Path c : children)
			c.setChildren(pathsLeft);
	}

	/**
	 * Returns the ratio between the "Euclidean distance" of this path and its
	 * length. The Euclidean distance of this path is defined as the distance
	 * between this Path's start and end point.
	 *
	 * @return the Contraction of this Path, or {@code NaN} if this Path has no length
	 */
	public double getContraction() {
		final double length = getLength();
		if (length == 0) return Double.NaN;
		final double eDistance = getNode(size() - 1).distanceTo(getNode(0));
		return eDistance / length;
	}

	/**
	 * Gets the fractal dimension of this path.
	 *
	 * @return the fractal dimension of this path or Double.NaN if this path has less than 5 nodes
	 */
	public double getFractalDimension() {
		if (size() < 5) { // Must have at least 4 points after the start-node
			return Double.NaN;
		}
		final List<Double> pathDists = new ArrayList<>();
		final List<Double> eucDists = new ArrayList<>();
		for (int i = 1; i < size(); i++) { //Start at the second node
			double pDist = getNodeWithoutChecks(i).distanceTo(getNodeWithoutChecks(i - 1));
			if (!pathDists.isEmpty()) {
				double cumDist = pathDists.get(i - 2) + pDist;
				pathDists.add(cumDist);
			} else {
				pathDists.add(pDist);
			}
			double eDist = getNodeWithoutChecks(i).distanceTo(getNodeWithoutChecks(0));
			eucDists.add(eDist);
		}
		double numerator = 0.0;
		for (int i = 0; i < eucDists.size(); i++) {
			numerator += Math.log(1 + eucDists.get(i)) * Math.log(1 + pathDists.get(i));
		}
		double denominator = 0.0;
		for (final double eucDist : eucDists) {
			denominator += Math.log(1 + eucDist) * Math.log(1 + eucDist);
		}
		return numerator / denominator;
	}

	/**
	 * Gets the length of this Path
	 *
	 * @return the length of this Path
	 */
	public double getLength() {
		if (nodes.size() < 2) return 0;
		
		double totalLength = 0;
		PathNode prevNode = nodes.get(0);
		for (int i = 1; i < nodes.size(); ++i) {
			final PathNode currentNode = nodes.get(i);
			final double xdiff = currentNode.x - prevNode.x;
			final double ydiff = currentNode.y - prevNode.y;
			final double zdiff = currentNode.z - prevNode.z;
			totalLength += Math.sqrt(xdiff * xdiff + ydiff * ydiff + zdiff * zdiff);
			prevNode = currentNode;
		}
		return totalLength;
	}

	/**
	 * Computes the angle between the specified node and its two flanking neighbors.
	 * <p>
	 * With B being the specified node, A its previous neighbor, and C is next neighbor, computes the angle
	 * between the vectors AB, and BC.
	 * </p>
	 * @return the angle in degrees (0-360 range) or Double.NaN if specified node does not have sufficient neighbors
	 */
	public double getAngle(final int nodeIndex) {
		if (nodeIndex < 3)
			return Double.NaN;
		final PointInImage a = getNode(nodeIndex - 1);
		final PointInImage b = getNode(nodeIndex - 2);
		final PointInImage c = getNode(nodeIndex);
		final Vector3d v1 = new Vector3d(a.getX() - b.getX(), a.getY() - b.getY(), a.getZ() - b.getZ()); // BA
		final Vector3d v2 = new Vector3d(c.getX() - b.getX(), c.getY() - b.getY(), c.getZ() - b.getZ()); // BC
		v1.normalize();
		v2.normalize();
		double angle = Math.acos(v1.dot(v2));
		angle = MathUtils.normalizeAngle(angle, Math.PI); // normalize angle between 0 and 2PI
		return Math.toDegrees(angle);
	}

	/**
	 * Returns the overall extension angle of this path in the XY plane.
	 * The angle is obtained from the slope of a linear regression across all the path nodes.
	 *
	 * @return the overall 'extension' angle in degrees [0-360[ of this path in the XY plane.
	 */
	public double getExtensionAngleXY() {
		return getExtensionAngle(MultiDThreePanes.XY_PLANE, false);
	}

	/**
	 * Returns the overall extension angle of this path in the XZ plane.
	 * The angle is obtained from the slope of a linear regression across all the path nodes.
	 *
	 * @return the overall 'extension' angle in degrees [0-360[ in the XZ plane or NaN if path is 2D.
	 */
	public double getExtensionAngleXZ() {
		return (is3D()) ? getExtensionAngle(MultiDThreePanes.XZ_PLANE, false) : Double.NaN;
	}

	/**
	 * Returns the overall extension angle of this path in the ZY plane.
	 * The angle is obtained from the slope of a linear regression across all the path nodes.
	 *
	 * @return the overall 'extension' angle in degrees [0-360[ in the ZY plane or NaN if path is 2D.
	 */
	public double getExtensionAngleZY() {
		return (is3D()) ? getExtensionAngle(MultiDThreePanes.ZY_PLANE, false) : Double.NaN;
	}

	private boolean is3D() throws IllegalArgumentException {
		final double zRef = getNodeWithoutChecks(0).getZ();
		for (int i=0; i< size(); i++) {
			if (getNodeWithoutChecks(i).getZ() != zRef) return true;
		}
		return false;
	}

	private double getExtensionAngle(final int view, final boolean relative) {
		if (!relative)
			return getExtensionAngle(view);
		if (parentPath == null)
			return Double.NaN;
		return getExtensionAngle3D(parentPath.getExtensionDirection3D());
	}

	private double getExtensionAngle(final int view) {
		if (size() == 1)
			return Double.NaN;
        double slope;
        if (size() == 2) {
            slope = switch (view) {
                case MultiDThreePanes.XY_PLANE -> getNodeWithoutChecks(1).y - getNodeWithoutChecks(0).y
                        / getNodeWithoutChecks(1).x - getNodeWithoutChecks(0).x;
                case MultiDThreePanes.XZ_PLANE -> getNodeWithoutChecks(1).z - getNodeWithoutChecks(0).z
                        / getNodeWithoutChecks(1).x - getNodeWithoutChecks(0).x;
                case MultiDThreePanes.ZY_PLANE -> getNodeWithoutChecks(1).y - getNodeWithoutChecks(0).y
                        / getNodeWithoutChecks(1).z - getNodeWithoutChecks(0).z;
                default -> throw new IllegalArgumentException("Not a valid plane");
            };
        } else {
            final SimpleRegression sr = new SimpleRegression();
            switch (view) {
                case MultiDThreePanes.XY_PLANE:
                    getNodes().forEach(node -> sr.addData(node.x, node.y));
                    break;
                case MultiDThreePanes.XZ_PLANE:
                    getNodes().forEach(node -> sr.addData(node.x, node.z));
                    break;
                case MultiDThreePanes.ZY_PLANE:
                    getNodes().forEach(node -> sr.addData(node.z, node.y));
                    break;
                default:
                    throw new IllegalArgumentException("Not a valid plane");
            }
            slope = sr.getSlope(); // Get the slope from the regression line
        }

		// Handle vertical lines (infinite slope or NaN slope) separately
		if (Double.isInfinite(slope) || Double.isNaN(slope)) {
			// For vertical lines, determine direction by looking at Y coordinates
			double y1, y2;
			switch (view) {
				case MultiDThreePanes.XY_PLANE, MultiDThreePanes.ZY_PLANE:
					y1 = getNodeWithoutChecks(0).y;
					y2 = getNodeWithoutChecks(size() - 1).y;
					break;
				case MultiDThreePanes.XZ_PLANE:
					y1 = getNodeWithoutChecks(0).z;
					y2 = getNodeWithoutChecks(size() - 1).z;
					break;
                default:
					y1 = y2 = 0;
			}
			// In image coordinates, Y increases downward, so North is negative Y direction
			if (y2 <= y1) {
				return 0.0; // North
			} else {
				return 180.0; // South
			}
		}
		
		// Calculate direction vector from regression line. We need to determine the direction
		// (positive or negative X direction) by looking at the overall X range of the data
		double x1, x2;
		switch (view) {
			case MultiDThreePanes.XY_PLANE, MultiDThreePanes.XZ_PLANE:
				x1 = getNodeWithoutChecks(0).x;
				x2 = getNodeWithoutChecks(size() - 1).x;
				break;
            case MultiDThreePanes.ZY_PLANE:
				x1 = getNodeWithoutChecks(0).z;
				x2 = getNodeWithoutChecks(size() - 1).z;
				break;
			default:
				x1 = x2 = 0; // This won't be reached due to earlier validation
		}
		
		// Determine direction: if path goes from left to right, deltaX = 1, else deltaX = -1
		final double deltaX = (x2 >= x1) ? 1.0 : -1.0;
		final double deltaY = slope * deltaX;
		double angle;

		// Calculate angle using atan2 for full quadrant information
		// Convert to navigation/compass convention: North: 0°; East: 90°; South: 180°; West: 270°
		// In image coordinates, Y increases downward, so North has a negative deltaY
		angle = Math.atan2(-deltaY, deltaX); // Flip Y axis for proper North direction
		
		// Now convert from standard math angles to compass bearing
		// Math: East=0°, North=90°; Compass: North=0°, East=90°
		angle = Math.PI / 2 - angle;
		
		// Normalize angle to [0, 2π] range
		angle = MathUtils.normalizeAngle(angle, Math.PI);
		return Math.toDegrees(angle); // return angle in 0-360 degrees
	}

	/**
	 * Returns a (non-normalized) direction vector from the first node to the
	 * node at index {@code min(maxNodes, size()-1)}. This is useful for
	 * measuring the initial heading of a child path near a fork point.
	 *
	 * @param maxNodes the maximum node index to use as the end point
	 * @return a 3-element array {@code [dx, dy, dz]}, or {@code null} if the
	 *         path has fewer than 2 nodes or the resulting vector is zero-length
	 */
	public double[] getInitialDirection(final int maxNodes) {
		final int end = Math.min(maxNodes, size() - 1);
		if (end < 1) return null;
		final PointInImage p0 = getNodeWithoutChecks(0);
		final PointInImage pN = getNodeWithoutChecks(end);
		final double dx = pN.x - p0.x;
		final double dy = pN.y - p0.y;
		final double dz = pN.z - p0.z;
		if (dx == 0 && dy == 0 && dz == 0) return null;
		return new double[]{ dx, dy, dz };
	}

	/**
	 * Returns the 3D extension direction vector of this path using linear regression.
	 * The vector represents the overall direction from start to end of the path, fitted through all its nodes.
	 *
	 * @return the normalized Vector3d (length of 1) representing the extension direction,
	 * or null if path has only one point
	 */
	public Vector3d getExtensionDirection3D() {
		if (size() == 1)
			return null;

        if (size() == 2) {
            final Vector3d v = getNodeWithoutChecks(1).asVector();
            v.sub(getNodeWithoutChecks(0).asVector());
            return v;
        }

		// Parametric approach: Use linear regression for each coordinate
		final SimpleRegression xRegression = new SimpleRegression();
		final SimpleRegression yRegression = new SimpleRegression();
		final SimpleRegression zRegression = new SimpleRegression();

		// Add data points using index as parameter
		for (int i = 0; i < size(); i++) {
			final PointInImage node = getNodeWithoutChecks(i);
			xRegression.addData(i, node.x);
			yRegression.addData(i, node.y);
			zRegression.addData(i, node.z);
		}

		// Calculate direction vector using regression slopes. Because we use indices
		// this is the same as xRegression.predict(size() - 1) - xRegression.predict(0)
		// but more direct and slightly more efficient:
		// given y = mx + b:
		//    predict(size() - 1) = m * (size() - 1) + b
		//    predict(0) = m * 0 + b = b
		//    predict(size() - 1) - predict(0) = m * (size() - 1) + b - b = m * (size() - 1)
		final double deltaX = xRegression.getSlope() * (size() - 1);
		final double deltaY = yRegression.getSlope() * (size() - 1);
		final double deltaZ = zRegression.getSlope() * (size() - 1);

		final Vector3d direction = new Vector3d(deltaX, deltaY, deltaZ);
		direction.normalize(); // Convert to unit vector
		return direction;
	}

	/**
	 * Returns either the compass bearing or the relative branching angle of this path's 3D extension direction.
	 * 
	 * @param relative determines the calculation mode:
	 *                 - false: returns compass bearing in XY plane (absolute angle)
	 *                 - true: returns angle relative to parent path (relative angle)
	 * @return when relative=false: compass bearing in degrees (0-360°) using navigation convention:
	 *         0° = North (negative Y direction in image coordinates),
	 *         90° = East (positive X direction),
	 *         180° = South (positive Y direction in image coordinates),
	 *         270° = West (negative X direction).
	 *         The 3D direction vector is projected onto the XY plane for compass calculation.
	 *         <p>
	 *         When relative=true: acute angle in degrees (0-180°) between this path's 3D direction
	 *         vector and its parent path's 3D direction vector. This measures the true 3D branching
	 *         angle, not just the XY projection.
	 *         <p>
	 *         Returns Double.NaN if:
	 *         - Path has only one point (cannot determine direction)
	 *         - relative=true and path has no parent (startJoins is null)
	 *         - relative=true and parent path has no direction (single point parent)
	 * @see #getExtensionDirection3D()
	 * @see #getExtensionAngle3D(Vector3d)
	 * @see #getExtensionAngles3D()
	 * @see #getExtensionAngleXY()
	 */
	public double getExtensionAngle3D(final boolean relative) {
		if (relative) {
			if (parentPath == null) return Double.NaN;
			final Vector3d sjDirection = parentPath.getExtensionDirection3D();
			return (sjDirection == null) ? Double.NaN : getExtensionAngle3D(sjDirection);
		}

		final Vector3d direction = getExtensionDirection3D();
		if (direction == null)
			return Double.NaN;
		// Use X and Y components for compass calculation
		final double deltaX = direction.x;
		final double deltaY = direction.y;

		// Calculate angle using atan2 for full quadrant information
		// Convert to navigation/compass convention: North: 0°; East: 90°; South: 180°; West: 270°
		// In image coordinates, Y increases downward, so North has a negative deltaY
		double angle = Math.atan2(-deltaY, deltaX); // Flip Y axis for proper North direction
		
		// Now convert from standard math angles to compass bearing
		// Math: East=0°, North=90°; Compass: North=0°, East=90°
		angle = Math.PI / 2 - angle;
		
		// Normalize angle to [0, 2π] range
		angle = MathUtils.normalizeAngle(angle, Math.PI);
		return Math.toDegrees(angle); // return angle in 0-360 degrees
	}

	/**
	 * Returns the complete 3D orientation of this path's extension direction as spherical coordinates
	 * using navigation/compass convention.
	 * <p>
	 * This method provides both horizontal direction (azimuth) and vertical inclination (elevation)
	 * of the path's overall extension direction.
	 * 
	 * @return double array containing [azimuth, elevation] in degrees, where:
	 *         <p>
	 *         <b>azimuth</b> (index 0): Compass bearing in XY plane (0-360°) following navigation convention:
	 *         <ul>
	 *         <li>0° = North (negative Y direction in image coordinates)</li>
	 *         <li>90° = East (positive X direction)</li>
	 *         <li>180° = South (positive Y direction in image coordinates)</li>
	 *         <li>270° = West (negative X direction)</li>
	 *         </ul>
	 *         <p>
	 *         <b>elevation</b> (index 1): Vertical angle from XY plane (-90° to +90°):
	 *         <ul>
	 *         <li>0° = horizontal (parallel to XY plane)</li>
	 *         <li>+90° = pointing straight up (positive Z direction)</li>
	 *         <li>-90° = pointing straight down (negative Z direction)</li>
	 *         </ul>
	 *         <p>
	 *         Returns null if path has only one point or extension direction cannot be determined
	 *         <p>
	 * 
	 * @see #getExtensionDirection3D()
	 * @see #getExtensionAngle3D(boolean)
	 * @see #getExtensionAngleXY()
	 * @see #getExtensionAngleFromVertical()
	 */
	public double[] getExtensionAngles3D() {
		final Vector3d direction = getExtensionDirection3D();
		if (direction == null)
			return null;
		
		// Calculate azimuth angle in XY plane using navigation/compass convention
		// 0°: North (negative Y), 90°: East (positive X), 180°: South (positive Y), 270°: West (negative X)
		// Note: In image coordinates, Y increases downward, so North is negative Y direction
		double azimuth = Math.atan2(direction.x, -direction.y);
		azimuth = MathUtils.normalizeAngle(azimuth, Math.PI); // normalize to [0, 2π]
		azimuth = Math.toDegrees(azimuth);
		
		// Calculate elevation angle from XY plane (-90° to +90°)
		final double xyLength = Math.sqrt(direction.x * direction.x + direction.y * direction.y);
		double elevation = Math.atan2(direction.z, xyLength);
		elevation = Math.toDegrees(elevation);
		
		return new double[]{azimuth, elevation};
	}

	/**
	 * Returns a single angle representing the 3D extension direction.
	 * This is the angle between the path's direction vector and a reference vector.
	 * 
	 * @param referenceVector the reference vector to measure angle from (e.g., new Vector3d(0, 0, 1) for vertical
	 *                       reference, or new Vector3d(1, 0, 0) for horizontal)
	 * @return the acute angle in degrees (0-180°), or Double.NaN if path has only one point
	 */
	public double getExtensionAngle3D(final Vector3d referenceVector) {
		final Vector3d direction = getExtensionDirection3D();
		if (direction == null)
			return Double.NaN;
		
		// Calculate angle between direction and reference vector
		final double dotProduct = direction.dot(referenceVector);
		// Clamp to [-1, 1] to handle floating point precision issues
		final double clampedDot = Math.max(-1.0, Math.min(1.0, dotProduct));
		final double angle = Math.acos(clampedDot);
		
		return Math.toDegrees(angle);
	}

	/**
	 * Returns the angle between the path's 3D direction and the vertical axis.
	 *
	 * @return angle from vertical in degrees (0° = vertical, 90° = horizontal), 
	 *         or Double.NaN if path has only one point, or 90° if path is 2D
	 */
	public double getExtensionAngleFromVertical() {
		return getExtensionAngle3D(new Vector3d(0, 0, 1)); // Use positive Z as vertical reference
	}

	/**
	 * Returns the angle between the path's 3D direction and the horizontal axis.
	 *
	 * @return angle from horizontal in degrees (09° = vertical, 0° = horizontal),
	 *         or Double.NaN if path has only one point.
	 */
	public double getExtensionAngleFromHorizontal() {
		return getExtensionAngle3D(new Vector3d(1, 0, 0)); // Use positive X as horizontal reference
	}

	protected String getRealLengthString() {
		return String.format(Locale.US, "%.3f", getLength()); // see https://github.com/morphonets/SNT/issues/147
	}

	/**
	 * @deprecated No longer needed - radii and tangents are now stored per-node
	 */
	@Deprecated
	public void createCircles() {
		// No-op in new design - radii and tangents are stored per-node as needed
	}

	protected void setIsPrimary(final boolean primary) {
		if (primary) setOrder(1);
	}

	/**
	 * Checks if this Path is root.
	 *
	 * @return true, if is primary (root)
	 */
	public boolean isPrimary() {
		return order == 1 || (parentPath == null);
	}

	/*
	 * We call this if we're going to delete the path represented by this object
	 */
	protected void disconnectFromAll() {
		/*
		 * This path can be connected to other ones either if: 1) this starts on other;
		 * 2) other starts on this In any of these cases, we need to also remove this
		 * from other's somehowJoins and other from this's somehowJoins.
		 */
		for (final Path other : connectedPaths) {
			if (other.parentPath == this) {
				other.parentPath = null;
				other.branchPointIndex = -1;
			}
			other.connectedPaths.remove(this);
			other.children.remove(this);
		}
		connectedPaths.clear();
		children.clear();
		parentPath = null;
		branchPointIndex = -1;
		setIsPrimary(true);
	}

	/**
	 * Establishes this path as a child path originating from the specified parent path.
	 * <p>
	 * In detail:
	 * <ul>
	 * <li>Sets the parent path that this path branches from</li>
	 * <li>Records the exact point where branching occurs</li>
	 * <li>Adds this path to the parent's children collection</li>
	 * <li>Updates this path's order to be parent order + 1</li>
	 * <li>Maintains bidirectional references for tree traversal</li>
	 * </ul>
	 *
	 * @param parentPath  the parent path that this path branches from. Must not be null.
	 *                    Some operations may expect it to be an existing path in the same
	 *                    tree structure
	 * @param branchPoint the exact 3D coordinate where this path branches from the parent.
	 *                    Must not be null and should correspond to a location on or near
	 *                    the parent path
	 * @throws IllegalArgumentException if parentPath is null
	 * @throws IllegalArgumentException if this path already has a parent (use
	 *                                  {@link #detachFromParent()} first to detach)
	 * @see #getParentPath() to retrieve the parent path
	 * @see #getBranchPoint() to get the branching point
	 * @see #getChildren() to access child paths branching from this path
	 * @see #getOrder() for the hierarchical level (parent order + 1)
	 * @see #detachFromParent() for disconnecting this path
	 */
	public void setBranchFrom(final Path parentPath, final PointInImage branchPoint) {
		if (parentPath == null) {
			throw new IllegalArgumentException("setBranchFrom should never take a null path");
		}
		{
			// If there was an existing path, that's an error:
			if (this.parentPath != null)
				throw new IllegalArgumentException("setBranchFrom should not replace existing child-parent connections");
			this.parentPath = parentPath;
			
			// Find the nearest node index in the parent path
			if (branchPoint != null) {
				this.branchPointIndex = parentPath.indexNearestTo(branchPoint.x, branchPoint.y, branchPoint.z);
				if (this.branchPointIndex == -1) {
					// If no nearby node found, use the last node as fallback
					this.branchPointIndex = parentPath.size() - 1;
				}
			} else {
				this.branchPointIndex = -1;
			}
		}
		// Also update the somehowJoins list:
		if (!connectedPaths.contains(parentPath)) {
			connectedPaths.add(parentPath);
		}
		if (!parentPath.connectedPaths.contains(this)) {
			parentPath.connectedPaths.add(this);
		}
		if (!parentPath.children.contains(this)) {
			parentPath.children.add(this);
		}
		// update order
		setOrder(parentPath.order + 1);
	}

	/**
	 * @deprecated use {@link #setBranchFrom(Path, PointInImage)} instead
	 */
	@Deprecated
	public void setStartJoin(final Path parentPath, final PointInImage branchPoint) {
		setBranchFrom(parentPath, branchPoint);
	}

    protected void replaceNodesWithFittedVersion() {
        if (fitted == null) {
            throw new IllegalArgumentException("Path has no fitted version");
        }
        replaceNodes(fitted);
        if (!getName().contains(" [Fitted*]")) {
            setName(getName() + " [Fitted*]");
        }
        setUseFitted(false);
        fitted = null;
    }

	protected void rebuildConnectionsOfFittedVersion() {
		if (fitted == null)
			return;
		if (getParentPath() != null) { // this is always the case if not primary
			if (fitted.parentPath != null) fitted.detachFromParent();
			final PointInImage currentBranchPoint = getBranchPoint();
			if (currentBranchPoint != null) {
				if (parentPath.getUseFitted()) {
					fitted.setBranchFrom(parentPath.getFitted(), currentBranchPoint);
				} else {
					fitted.setBranchFrom(parentPath, currentBranchPoint);
				}
			}
		}
	}

	@Deprecated
	public void unsetStartJoin() {
		detachFromParent();
	}

	/**
	 * Detaches this path from its parent, converting it into an independent primary path.
	 * <p>
	 * Removes the parent-child relationship established by {@link #setBranchFrom(Path, PointInImage)}
	 * In
	 * <ul>
	 * <li>Clears the parent path reference and branch point</li>
	 * <li>Removes this path from the parent's children collection</li>
	 * <li>Updates bidirectional references for tree traversal</li>
	 * <li>Resets the path order to -1 (no hierarchical position)</li>
	 * <li>Converts this path from a branch to a primary path</li>
	 * </ul>
	 *
	 * @throws IllegalArgumentException if this path has no parent (i.e., is already
	 *                                  a primary path). Use {@link #getParentPath()}
	 *                                  to check before calling this method.
	 *
	 * @see #setBranchFrom(Path, PointInImage) to establish a parent-child relationship
	 * @see #getParentPath() to check if this path has a parent
	 * @see #getBranchPoint() to get the current branch point (before detaching)
	 * @see #isPrimary() to check if this path is already independent
	 * @see #getChildren() to access paths that branch from this path
	 * @see #getOrder() for hierarchical position information
	 */
	public void detachFromParent() {
		Path other = parentPath;
		if (other == null) {
			throw new IllegalArgumentException("Don't call detachFromParent if the other Path is already null");
		}
		if (!(other.parentPath == this)) {
			connectedPaths.remove(other);
			other.connectedPaths.remove(this);
			other.children.remove(this);
		}
		parentPath = null;
		branchPointIndex = -1;
		setOrder(-1);
	}

	public double getMinimumSeparation() {
		return Math.min(Math.abs(x_spacing), Math.min(Math.abs(y_spacing), Math.abs(z_spacing)));
	}

	/**
	 * Returns the number of nodes of this path
	 *
	 * @return the size, i.e., number of nodes
	 */
	public int size() {
		return nodes.size();
	}

	/**
	 * Gets all nodes as a list.
	 *
	 * @return a new list containing all nodes
	 */
	public List<PointInImage> getNodes() {
		return new ArrayList<>(nodes);
	}

	protected void getPointDouble(final int i, final double[] p) {
		if ((i < 0) || i >= size()) {
			throw new IndexOutOfBoundsException(
				"BUG: getPointDouble was asked for an out-of-range point: " + i);
		}

		final PathNode node = getNodeWithoutChecks(i);
		p[0] = node.x;
		p[1] = node.y;
		p[2] = node.z;
	}

	/**
	 * Gets the node at the specified position.
	 *
	 * @param pos the position index (0-based). NB: You can use -1 to retrieve the
	 *            last node in the Path
	 * @return the node
	 * @throws IndexOutOfBoundsException if position is out-of-range
	 */
	public PathNode getNode(final int pos) throws IndexOutOfBoundsException {
		return getNodeWithChecks(pos);
	}

    protected PathNode getNode(String filter) {
        switch (filter.toLowerCase().trim()) {
            case "min radius" -> {
                if (!hasRadii()) return null;
                return nodes.stream()
                        .min(java.util.Comparator.comparingDouble(Path.PathNode::getRadius))
                        .orElse(null);
            }
            case "max radius" -> {
                if (!hasRadii()) return null;
                return nodes.stream()
                        .max(java.util.Comparator.comparingDouble(Path.PathNode::getRadius))
                        .orElse(null);
            }
            case "min x" -> {
                return nodes.stream()
                        .min(java.util.Comparator.comparingDouble(Path.PathNode::getX))
                        .orElse(null);
            }
            case "max x" -> {
                return nodes.stream()
                        .max(java.util.Comparator.comparingDouble(Path.PathNode::getX))
                        .orElse(null);
            }
            case "min y" -> {
                return nodes.stream()
                        .min(java.util.Comparator.comparingDouble(Path.PathNode::getY))
                        .orElse(null);
            }
            case "max y" -> {
                return nodes.stream()
                        .max(java.util.Comparator.comparingDouble(Path.PathNode::getY))
                        .orElse(null);
            }
            case "min z" -> {
                return nodes.stream()
                        .min(java.util.Comparator.comparingDouble(Path.PathNode::getZ))
                        .orElse(null);
            }
            case "max z" -> {
                return nodes.stream()
                        .max(java.util.Comparator.comparingDouble(Path.PathNode::getZ))
                        .orElse(null);
            }
            default -> {
                throw new IllegalArgumentException("unsupported filter");
            }
        }
    }

	/**
	 * Gets the internal PathNode at the specified position (for advanced operations).
	 */
	protected PathNode getNodeWithoutChecks(final int pos) {
		return nodes.get(pos);
	}

	/**
	 * Gets the internal PathNode at the specified position (for advanced operations).
	 */
	protected PathNode getNodeWithChecks(final int pos) {
		if (pos == -1)
			return getNodeWithoutChecks(size()-1);
		if ((pos < 0) || pos >= size()) {
			throw new IndexOutOfBoundsException("getNode() was asked for an out-of-range point: " + pos);
		}
		return nodes.get(pos);
	}

	public PointInCanvas getPointInCanvas(final int node) {
		final PointInCanvas result = new PointInCanvas(getXUnscaledDouble(node), getYUnscaledDouble(node),
				getZUnscaledDouble(node));
		result.onPath = this;
		return result;
	}

	protected PointInCanvas getPointInCanvas2D(final int node, final TracerCanvas canvas) {
		final PointInCanvas result = switch (canvas.getPlane()) {
            case MultiDThreePanes.XZ_PLANE -> new PointInCanvas(getXUnscaledDouble(node), getZUnscaledDouble(node), 0);
            case MultiDThreePanes.ZY_PLANE -> new PointInCanvas(getZUnscaledDouble(node), getYUnscaledDouble(node), 0);
            default -> new PointInCanvas(getXUnscaledDouble(node), getYUnscaledDouble(node), 0);
        };
        result.onPath = this;
		return result;
	}

	protected List<PointInCanvas> getUnscaledNodesInViewPort(final TracerCanvas canvas) {
		final BoundingBox box = canvas.getViewPort();
		final List<PointInCanvas> nodes = new ArrayList<>();
		for (int i = 0; i < size(); i++) {
			final PointInCanvas node = getPointInCanvas2D(i, canvas);
			if (box.contains(node)) nodes.add(node);
		}
		return nodes;
	}

	protected boolean containsUnscaledNodesInViewPort(final TracerCanvas canvas) {
		final BoundingBox box = canvas.getViewPort();
		for (int i = 0; i < size(); i++) {
			final PointInCanvas node = getPointInCanvas2D(i, canvas);
			if (box.contains(node))
				return true;
		}
		return false;
	}

	/**
	 * Check whether this Path contains the specified point
	 *
	 * @param pim the {@link PointInImage} node
	 * @return true, if successful
	 */
	public boolean contains(final PointInImage pim) {
		if (pim.onPath != null) return pim.onPath == this;
		// Direct iteration is faster than stream for simple checks
		for (PathNode node : nodes) {
			if (node.x == pim.x && node.y == pim.y && node.z == pim.z) {
				return true;
			}
		}
		return false;
	}

	/**
	 * Inserts a node at a specified position.
	 *
	 * @param index the (zero-based) index of the position of the new node
	 * @param point the node to be inserted
	 * @throws IllegalArgumentException if index is out-of-range
	 */
	public void insertNode(final int index, final PointInImage point) {
		if (index < 0 || index > size()) throw new IllegalArgumentException(
			"insertNode() asked for an out-of-range point: " + index);

		PathNode node = new PathNode(point);
		node.onPath = this;
		nodes.add(index, node);
	}

	/**
	 * Inserts a node at the specified position.
	 *
	 * @param index the position at which the node is to be inserted
	 * @param point the node to be inserted
	 * @throws IllegalArgumentException if index is out-of-range
	 */
	public void insertNode(final int index, final SNTPoint point) {
		insertNode(index, new PointInImage(point.getX(), point.getY(), point.getZ()));
	}

	/**
	 * Removes the specified node if this path has at least two nodes. Does
	 * nothing if this is a single point path.
	 *
	 * @param index the zero-based index of the node to be removed
	 * @throws IllegalArgumentException if index is out-of-range
	 */
	public void removeNode(final int index) {
		if (size() == 1) return;
		nodes.remove(index);
		
		// Update branch point indices for child paths that may be affected
		for (final Path child : getChildren()) {
			if (child.branchPointIndex > index) {
				child.branchPointIndex--; // Shift index down due to removal
			} else if (child.branchPointIndex == index) {
				// The branch point node was removed, find nearest remaining node
				if (!nodes.isEmpty()) {
					final PointInImage firstChildNode = child.firstNode();
					if (firstChildNode != null) {
						child.branchPointIndex = indexNearestTo(firstChildNode.x, firstChildNode.y, firstChildNode.z);
						if (child.branchPointIndex == -1) {
							child.branchPointIndex = Math.min(index, size() - 1);
						}
					}
				}
			}
		}
	}

	/**
	 * Assigns a new location to the specified node.
	 *
	 * @param index the zero-based index of the node to be modified
	 * @param destination the new node location
	 * @throws IllegalArgumentException if index is out-of-range
	 */
	public void moveNode(final int index, final PointInImage destination) {
		final PathNode node = getNodeWithChecks(index);
        // Update node coordinates
		node.x = destination.x;
		node.y = destination.y;
		node.z = destination.z;
		
		// Update fitted path if it exists (avoid recursion by checking if this is not already a fitted path)
		if (getFitted() != null && !isFittedVersionOfAnotherPath() && index < getFitted().size()) {
			getFitted().moveNode(index, destination);
		}
		
		// Note: Branch points now automatically stay in sync since they reference the actual parent nodes
		// No manual synchronization needed!
	}

    /**
     * Assigns a new location to the specified node using raw coordinates.
     * This is a performance-optimized version of {@link #moveNode(int, PointInImage)}
     * that avoids object allocation.
     *
     * @param index the index of the node to move
     * @param x the new x coordinate
     * @param y the new y coordinate
     * @param z the new z coordinate
     * @throws IndexOutOfBoundsException if the index is out of range
     */
    public void moveNode(final int index, final double x, final double y, final double z) {
        final PathNode node = getNodeWithChecks(index);
        node.x = x;
        node.y = y;
        node.z = z;
        if (getFitted() != null && !isFittedVersionOfAnotherPath() && index < getFitted().size()) {
            getFitted().moveNode(index, x, y, z);
        }
    }

	/**
	 * Gets the first node index associated with the specified image coordinates.
	 * Returns -1 if no such node was found.
	 *
	 * @param pim the image position (calibrated coordinates)
	 * @return the index of the first node occurrence or -1 if there is no such
	 *         occurrence
	 */
	public int getNodeIndex(final PointInImage pim) {
		// Direct list access is faster than method calls
		for (int i = 0; i < nodes.size(); ++i) {
			final PathNode node = nodes.get(i);
			if (Math.abs(node.x - pim.x) < x_spacing && 
				Math.abs(node.y - pim.y) < y_spacing && 
				Math.abs(node.z - pim.z) < z_spacing) {
				return i;
			}
		}
		return -1;
	}

	/**
	 * Gets the closest node associated with the specified node.
	 *
	 * @param node   location to be searched for
	 * @param within sets the search sensitivity. E.g., Setting it to
	 *               Double.MAX_VALUE (or the image's largest dimension) will always
	 *               return a valid node.
	 * @return the closest node to the specified one or null if no such node was
	 *         found.
	 */
	public PointInImage nearestNodeTo(final PointInImage node, final double within) {
		double minimumDistanceSquared = within * within;
		PointInImage closestNode = null;
		for (int i = 0; i < size(); ++i) {
			final PointInImage cNode = getNodeWithoutChecks(i);
			final double thisDistanceSquared = cNode.distanceSquaredTo(node);
			if (thisDistanceSquared < minimumDistanceSquared) {
				closestNode = cNode;
				minimumDistanceSquared = thisDistanceSquared;
			}
		}
		return closestNode;
	}

	/**
	 * Gets the index of the closest node associated with the specified world
	 * coordinates.
	 *
	 * @param x the x-coordinates (spatially calibrated image units)
	 * @param y the y-coordinates (spatially calibrated image units)
	 * @param z the z-coordinates (spatially calibrated image units)
	 * @param within sets the search sensitivity. E.g., Setting it to
	 *          Double.MAX_VALUE (or the image's largest dimension) will always
	 *          return a valid index.
	 * @return the index of the closest node to the specified coordinates. Returns
	 *         -1 if no such node was found.
	 */
	public int indexNearestTo(final double x, final double y, final double z,
		final double within)
	{

		if (nodes.isEmpty()) throw new IllegalArgumentException(
			"indexNearestTo called on a Path of size() = 0");

		double minimumDistanceSquared = within * within;
		int indexOfMinimum = -1;

		// Direct list access is faster than method calls
		for (int i = 0; i < nodes.size(); ++i) {
			final PathNode node = nodes.get(i);
			final double diff_x = x - node.x;
			final double diff_y = y - node.y;
			final double diff_z = z - node.z;

			final double thisDistanceSquared = diff_x * diff_x + diff_y * diff_y +
				diff_z * diff_z;

			if (thisDistanceSquared < minimumDistanceSquared) {
				indexOfMinimum = i;
				minimumDistanceSquared = thisDistanceSquared;
			}
		}

		return indexOfMinimum;
	}

	protected int indexNearestToCanvasPosition2D(final double xInCanvas,
		final double yInCanvas, final double withinPixels)
	{
		double minimumDistanceSquared = withinPixels * withinPixels;
		int indexOfMinimum = -1;
		for (int i = 0; i < size(); ++i) {
			final double diff_x = xInCanvas - getXUnscaledDouble(i);
			final double diff_y = yInCanvas - getYUnscaledDouble(i);
			final double thisDistanceSquared = diff_x * diff_x + diff_y * diff_y;
			if (thisDistanceSquared < minimumDistanceSquared) {
				indexOfMinimum = i;
				minimumDistanceSquared = thisDistanceSquared;
			}
		}
		return indexOfMinimum;
	}

	protected int indexNearestTo(final double x, final double y, final double z) {
		return indexNearestTo(x, y, z, Double.MAX_VALUE);
	}

	/**
	 * Gets the position of the node tagged as 'editable', if any.
	 *
	 * @return the index of the point currently tagged as editable, or -1 if no
	 *         such point exists
	 */
	public int getEditableNodeIndex() {
		return editableNodeIndex;
	}

	/**
	 * Tags the specified point position as 'editable'.
	 *
	 * @param index the index of the point to be tagged. Set it to -1 to for no
	 *          tagging
	 */
	public void setEditableNode(final int index) {
		this.editableNodeIndex = index;
	}

	protected boolean isBeingEdited() {
		return editableNodeIndex != -1;
	}

	protected boolean isEditableNodeLocked() {
		return editableNodeLocked;
	}

	protected void setEditableNodeLocked(final boolean editableNodeLocked) {
		this.editableNodeLocked = editableNodeLocked;
	}

	/**
	 * Returns the unscaled X position of a node.
	 *
	 * @param i the index of the path node
	 * @return the 0-based index of the X-position of the node.
	 */
	public int getXUnscaled(final int i) {
		return (int) Math.round(getXUnscaledDouble(i));
	}

	/**
	 * Returns the unscaled Y position of a node.
	 *
	 * @param i the index of the path node
	 * @return the 0-based index of the Y-position of the node.
	 */
	public int getYUnscaled(final int i) {
		return (int) Math.round(getYUnscaledDouble(i));
	}

	/**
	 * Returns the unscaled Z position of a node.
	 *
	 * @param i the index of the path node
	 * @return the 0-based index of the Z-position of the node. Note that Z-indices in ImagePlus are 1-based,
	 * i.e., outside SNT the proper Z-slice for this position in an ImagePlus is {@code getZUnscaled(i) + 1}
	 */
	public int getZUnscaled(final int i) { return (int) Math.round(getZUnscaledDouble(i)); }

	public double getXUnscaledDouble(final int i) {
		return nodes.get(i).x / x_spacing + canvasOffset.x;
	}

	public double getYUnscaledDouble(final int i) {
		return nodes.get(i).y / y_spacing + canvasOffset.y;
	}

	public double getZUnscaledDouble(final int i) {
		return nodes.get(i).z / z_spacing + canvasOffset.z;
	}

    /**
     * @return the list of nodes in image space (pixel coordinates). This takes already takes into account
     * {@link #getCanvasOffset()}.
     */
    public List<PointInCanvas> getUnscaledNodes() {
        List<PointInCanvas> pics = new ArrayList<>(nodes.size());
        nodes.forEach(node -> {
            pics.add(new PointInCanvas(
                    node.x / x_spacing + canvasOffset.x,
                    node.y / y_spacing + canvasOffset.y,
                    node.z / z_spacing + canvasOffset.z));
        });
        return pics;
    }

	/**
	 * Gets a section between two nodes of this Path.
	 *
	 * @param startIndex the node index defining the first node of the section
	 * @param endIndex   the node index defining the last node of the section
	 * @return the section. Note that the sectioned Path will share of all of this
	 *         Path's properties but will not contain any information on junction
	 *         points or connectivity to other Paths.
	 */
	public Path getSection(final int startIndex, final int endIndex) {
		if (startIndex < 0 || endIndex > size() || endIndex < startIndex) {
			throw new IllegalArgumentException("Indices out of range!");
		}
		final Calibration cal = getCalibration();
		final Path sub = new Path(cal.pixelWidth, cal.pixelHeight, cal.pixelDepth, cal.getUnit());
		
		// Copy nodes from the specified range
		for (int i = startIndex; i <= endIndex; i++) {
			final PathNode originalNode = getNodeWithoutChecks(i);
			final PathNode newNode = originalNode.clone();
			newNode.onPath = sub;
			sub.nodes.add(newNode);
		}
		
		if (getFitted() != null)
			sub.setFitted(getFitted().getSection(startIndex, endIndex));
		applyCommonProperties(sub);
		return sub;
	}

	private void applyCommonProperties(final Path other) {
		other.setOrder(getOrder());
		other.setIsPrimary(isPrimary());
		other.setSWCType(getSWCType());
		other.setCTposition(getChannel(), getFrame());
		other.setEditableNode((getEditableNodeIndex() < other.size()) ? getEditableNodeIndex() : -1);
		other.setCanvasOffset(getCanvasOffset());
		other.setColor(getColor());
		other.setSpineOrVaricosityCount(getSpineOrVaricosityCount());
		// Directly set ID without triggering change events during cloning
		other.id = id;
		other.editableNodeIndex = editableNodeIndex;
		other.editableNodeLocked = editableNodeLocked;
	}

	/**
	 * Returns a new Path with this Path's attributes (e.g. spatial scale), but no nodes.
	 *
	 * @return the empty path
	 */
	public Path createPath() {
		final Calibration cal = getCalibration();
		final Path dup = new Path(cal.pixelWidth, cal.pixelHeight, cal.pixelDepth, cal.getUnit());
		applyCommonProperties(dup);
		return dup;
	}

    /**
     * Replaces all nodes of this path with those from the replacement path,
     * preserving hierarchical relationships (parent, children, connections).
     * <p>
     * This method is useful for refining path geometry (e.g., from A* tracing)
     * while maintaining the tree structure.
     * </p>
     *
     * @param replacementPath the path providing the new node geometry
     * @throws IllegalArgumentException if replacementPath is null or empty
     */
    public void replaceNodes(final Path replacementPath) {
        if (replacementPath == null || replacementPath.size() == 0) {
            throw new IllegalArgumentException("Replacement path cannot be null or empty");
        }

        // Clear current nodes and copy from replacement path
        nodes.clear();
        for (int i = 0; i < replacementPath.size(); i++) {
            final PathNode replacementNode = replacementPath.getNodeWithoutChecks(i);
            final PathNode newNode = replacementNode.clone();
            newNode.onPath = this;
            nodes.add(newNode);
        }

        // Update branch point index to nearest node in parent (if has parent)
        if (parentPath != null && !nodes.isEmpty()) {
            final PathNode firstNode = nodes.getFirst();
            final int index = parentPath.indexNearestTo(firstNode.x, firstNode.y, firstNode.z);
            if (index != -1) {
                branchPointIndex = index;
            }
        }

        // Update children's branch point indices to reflect new geometry
        for (final Path child : children) {
            if (child.branchPointIndex >= 0) {
                final PointInImage childBranchPoint = child.getBranchPoint();
                if (childBranchPoint != null) {
                    final int newIndex = indexNearestTo(
                            childBranchPoint.x, childBranchPoint.y, childBranchPoint.z);
                    if (newIndex != -1) {
                        child.branchPointIndex = newIndex;
                    }
                }
            }
        }
    }

	protected void setSpacing(final Calibration cal) {
		this.x_spacing = cal.pixelWidth;
		this.y_spacing = cal.pixelHeight;
		this.z_spacing = cal.pixelDepth;
		this.spacing_units = cal.getUnit();
	}

    public void setSpacing(double xSpacing, double ySpacing, double zSpacing, String unit) {
        this.x_spacing = xSpacing;
        this.y_spacing = ySpacing;
        this.z_spacing = zSpacing;
        this.spacing_units = SNTUtils.getSanitizedUnit(unit);
    }

	/**
	 * Creates a deep copy of this Path.
	 * <p>
	 * This method creates a complete copy of the path including all nodes, properties,
	 * and calibration settings. The cloned path will have the same spatial coordinates,
	 * node properties (radius, color, tangent), and metadata as the original, but will
	 * be completely independent. Parent-child relationships are not preserved in the clone
	 * to avoid infinite recursion and should be established separately if needed.
	 * </p>
	 * <p>
	 * Fitted paths and legacy 3D viewer content are not cloned - the new path will need to
	 * be re-fitted and re-added to viewers if needed.
	 * </p>
	 *
	 * @return a new Path that is a deep copy of this path
	 * @see #clone(boolean) for the deprecated version with child inclusion option
	 */
	@SuppressWarnings("unchecked")
	@Override
	public Path clone() {
		// For Path, we use manual cloning due to final fields referencing mutable objects
		// Using super.clone() causes issues because the final 'nodes' field would reference
		// the same ArrayList instance, causing dup.nodes.clear() to clear the original's nodes
		
		final Calibration cal = getCalibration();
		final Path dup = new Path(cal.pixelWidth, cal.pixelHeight, cal.pixelDepth, cal.getUnit());

		// Clone all nodes
		for (final PathNode node : nodes) {
			PathNode clonedNode = node.clone();
			clonedNode.onPath = dup;
			dup.nodes.add(clonedNode);
		}

		// Clone path relationships (shallow clone of lists, but references remain)
		dup.connectedPaths = (ArrayList<Path>) connectedPaths.clone();
		dup.children = (ArrayList<Path>) children.clone();

		// Note: We don't deep clone parentPath, branchPointIndex, or fitted to avoid infinite recursion
		// These relationships should be established separately after cloning
		dup.parentPath = null;
		dup.branchPointIndex = -1;
		dup.fitted = null;

		// Clone canvas offset
		if (canvasOffset != null) {
			dup.canvasOffset = new PointInCanvas(canvasOffset.x, canvasOffset.y, canvasOffset.z);
			dup.canvasOffset.onPath = dup;
		}

		// Reset 3D viewer content (should not be shared)
		dup.content3D = null;
		dup.nameWhenAddedToViewer = null;
		
		// Clear change listeners to prevent events during cloning
		dup.changeListeners.clear();

		applyCommonProperties(dup);
		return dup;
	}

	/**
	 * Creates a copy of this Path with optional inclusion of immediate children.
	 * <p>
	 * This method creates a clone of this path and optionally includes clones of its
	 * immediate child paths. When children are included, they are properly reconnected
	 * to the cloned parent path, maintaining the parent-child relationships. Note that
	 * only immediate children are cloned - grandchildren and deeper descendants are not
	 * included.
	 * </p>
	 * <p>
	 * <strong>Limitations:</strong> This method only clones the immediate children and
	 * does not recursively clone the entire subtree. For complex tree structures, consider
	 * using {@link Tree#clone()} instead.
	 * </p>
	 *
	 * @param includeImmediateChildren if true, includes clones of immediate child paths;
	 *                                 if false, returns a clone without children
	 * @return a new Path that is a copy of this path, optionally including children
	 */
	public Path clone(final boolean includeImmediateChildren) {
		final Path dup = clone();
		if (!includeImmediateChildren) return dup;

		// Clone immediate children and reconnect them to the cloned parent
		dup.children = new ArrayList<>(children.size());
		for (final Path child : children) {
			final Path dupChild = child.clone();

			// Clear the cloned child's relationships since they point to original paths
			dupChild.connectedPaths.clear();
			dupChild.children.clear();

			// Reconnect the child to the cloned parent if it had a branch point
			if (child.getBranchPoint() != null) {
				dupChild.setBranchFrom(dup, child.getBranchPoint());
			}

			dup.children.add(dupChild);
		}
		return dup;
	}

	/**
	 * Gets the list of child paths that branch from this path.
	 * <p>
	 * Returns the collection of paths that have this path as their parent.
	 * The returned list is the actual internal list, so modifications to it
	 * will affect the path's structure.
	 * </p>
	 * 
	 * @return the List of child paths. May be empty if this path has no children,
	 *         but never null.
	 * 
	 * @see #setBranchFrom(Path, PointInImage) for creating parent-child relationships
	 * @see #getParentPath() for getting the parent of this path
	 * @see #isPrimary() for checking if this is a root path
	 */
	public List<Path> getChildren() {
		return children;
	}

	/**
	 * Gets the spatial calibration of this Path.
	 *
	 * @return the calibration details associated with this Path
	 */
	public Calibration getCalibration() {
		final Calibration cal = new Calibration();
		cal.setUnit(spacing_units);
		cal.pixelWidth = x_spacing;
		cal.pixelHeight = y_spacing;
		cal.pixelDepth = z_spacing;
		return cal;
	}

	public void add(final Path other) {

		if (other == null) {
			SNTUtils.warn("BUG: Trying to add null Path");
			return;
		}

		// If we're trying to add a path with radii to one that has none,
		// set default radii for existing nodes
		if (other.hasRadii() && !hasRadii()) {
			final double defaultRadius = getMinimumSeparation() * 2;
			for (final PathNode node : nodes)
				node.radius = defaultRadius;
		}

		int toSkip = 0;

		// Skip points at the beginning of the other path if they're
		// the same as the last point on this path
		if (!nodes.isEmpty() && !other.nodes.isEmpty()) {
			final PathNode lastNode = nodes.getLast();
			while (toSkip < other.size()) {
				final PathNode otherNode = other.nodes.get(toSkip);
				if (otherNode.x == lastNode.x && 
					otherNode.y == lastNode.y && 
					otherNode.z == lastNode.z) {
					toSkip++;
				} else {
					break;
				}
			}
		}

		// Add nodes from the other path, starting from toSkip
		for (int i = toSkip; i < other.size(); i++) {
			final PathNode otherNode = other.getNodeWithoutChecks(i);
			final PathNode newNode = otherNode.clone();
			newNode.onPath = this;
			nodes.add(newNode);
		}
		if (hasRadii()) setGuessedTangents(2);
	}

	/**
	 * Reverses this path so that its starting node becomes the last and vice versa.
	 */
	public void reverse() {
		reverseInSitu(this);
	}

	/**
	 * 
	 * @return a reversed version of this path in which node coordinates are
	 *         reversed. Note that for legacy reasons only the node coordinates are
	 *         reversed. Other properties (node colors, etc.) are not included in
	 *         the returned path
	 */
	public Path reversed() {
		final Path c = createPath();
		// Add nodes in reverse order
		for (int i = size() - 1; i >= 0; i--) {
			final PathNode originalNode = getNodeWithoutChecks(i);
			final PathNode reversedNode = originalNode.clone();
			reversedNode.onPath = c;
			c.nodes.add(reversedNode);
		}
		if (c.fitted != null) c.fitted = c.fitted.reversed();
		return c;
	}

	private static void reverseInSitu(final Path c) {
		Collections.reverse(c.nodes);
		if (c.editableNodeIndex > -1)
			c.editableNodeIndex = c.size() - 1 - c.editableNodeIndex;
		if (c.fitted != null)
			reverseInSitu(c.fitted);
	}

	/**
	 * Appends a node to this Path.
	 *
	 * @param point the node to be inserted
	 */
	public void addNode(final PointInImage point) {
		final PathNode node = new PathNode(point);
		node.onPath = this;
		nodes.add(node);
	}

	public void addPointDouble(final double x, final double y, final double z) {
		addNode(new PointInImage(x, y, z));
	}

	public void drawPathAsPoints(final TracerCanvas canvas, final Graphics2D g,
		final java.awt.Color c, final int plane, final boolean highContrast,
		final boolean drawDiameter)
	{
		drawPathAsPoints(canvas, g, c, highContrast, drawDiameter, 0, -1);
	}

	protected void drawPathAsPoints(final TracerCanvas canvas, final Graphics2D g,
		final java.awt.Color c, final int plane, final boolean drawDiameter,
		final int slice, final int either_side)
	{
		drawPathAsPoints(canvas, g, c, false, drawDiameter, slice,
			either_side);
	}

	protected void drawPathAsPoints(final Graphics2D g2,
		final TracerCanvas canvas, final SNT snt)
	{
		final boolean customColor = (hasCustomColor && snt.displayCustomPathColors);
		Color color = snt.deselectedColor;
		if (isSelected() && !customColor)
			color = snt.selectedColor;
		else if (customColor)
			color = getColor();
		final int sliceZeroIndexed = canvas.getImage().getZ() - 1;
		int eitherSideParameter = canvas.eitherSide;
		if (!canvas.just_near_slices)
			eitherSideParameter = -1;
		drawPathAsPoints(canvas, g2, color, customColor,
			snt.getDrawDiameters(), sliceZeroIndexed, eitherSideParameter);
	}

    public void drawPathAsPoints(final TracerCanvas canvas, final Graphics2D g2,
                                 final java.awt.Color c, final boolean highContrast,
                                 boolean drawDiameter, final int slice, final int either_side)
    {
        final int nodeCount = nodes.size();
        if (nodeCount == 0) return;

        // Pre-compute values constant across the loop
        final BasicStroke connectionStroke = new BasicStroke(
                (float) (canvas.nodeDiameter() / 3), BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND);
        final PointInImage branchPoint = getBranchPoint();
        final int editableNodeIndex = getEditableNodeIndex();
        final int defaultTransparency = canvas.getDefaultTransparency();
        final int outOfBoundsTransparency = canvas.getOutOfBoundsTransparency();

        int startIndexOfLastDrawnLine = -1;

        PathNodeCanvas previousNode = null;
        PathNodeCanvas currentNode = null;

        for (int i = 0; i < nodeCount; ++i) {

            previousNode = currentNode;  // Reuse from previous iteration
            currentNode = new PathNodeCanvas(this, i, canvas);

            // Only create a new previousNode for special cases
            if (i == 0 && branchPoint != null) {
                previousNode = new PathNodeCanvas(branchPoint, i, canvas);
            }

            final boolean outOfDepthBounds = (either_side >= 0) && (Math.abs(
                    currentNode.getSlice() - slice) > either_side);

            // Draw node
            if (!outOfDepthBounds) {
                currentNode.setEditable(editableNodeIndex == i);
                currentNode.draw(g2, c);
            }

            // Options for drawing inter-node segments and node diameter. Color
            // will be whatever has been set by PathNodeCanvas#draw(). If 2D canvas,
            // out-of-bounds transparency is ignored
            g2.setStroke(connectionStroke);
            g2.setColor(SNTColor.alphaColor(g2.getColor(),
                    (outOfDepthBounds) ? outOfBoundsTransparency : defaultTransparency));

            // We are within Z-bounds and have been asked to draw the diameters, just do it in XY
			if (drawDiameter && !outOfDepthBounds && !currentNode.skipDiameterForDirectionArrows()) {
				currentNode.drawDiameter(g2, slice, either_side);
			}

            // Draw connection to previous node
			if (previousNode != null) {
				if (startIndexOfLastDrawnLine != i - 1) {
					currentNode.drawConnection(g2, previousNode);
					startIndexOfLastDrawnLine = i - 1;
				}
				// Draw direction arrow (previous points toward current)
				previousNode.drawDirectionArrow(g2, currentNode);
			}

            // Handle branch point special case (first node)
            if (i == 0 && branchPoint != null) {
                final PathNodeCanvas jointNode = new PathNodeCanvas(branchPoint, 0, canvas);
                jointNode.setType(PathNodeCanvas.JOINT);
                jointNode.draw(g2, c);
                currentNode.setType(PathNodeCanvas.SLAB);
				currentNode.drawConnection(g2, previousNode);
            }
        }
    }

	/**
	 * Sets the node colors.
	 *
	 * @param colors the colors used to render the nodes of this. If null (the
	 *          default) all nodes are rendered using the Path color.
	 */
	public void setNodeColors(final Color[] colors) {
		if (colors == null) {
			nodes.forEach(node -> node.setColor(null));
		} else if (colors.length != size()) {
			throw new IllegalArgumentException("colors array must have as many elements as nodes");
		} else {
			for (int i = 0; i < colors.length; i++)
				getNodeWithoutChecks(i).setColor(colors[i]);
		}
	}

	/**
	 * Gets the node colors.
	 *
	 * @return the colors used to render the nodes of this path, or null if nodes
	 *         are rendered using the Path color
	 */
	public Color[] getNodeColors() {
		// Direct array creation is faster than using streams!?
		final Color[] colors = new Color[nodes.size()];
		for (int i = 0; i < nodes.size(); i++)
			colors[i] = nodes.get(i).getColor();
		return colors;
	}

	/**
	 * Gets the node color.
	 *
	 * @param pos the node position
	 * @return the node color, or null if no color nodes have been assigned to
	 *         this path
	 */
	public Color getNodeColor(final int pos) {
		return getNodeWithChecks(pos).getColor();
	}

	/**
	 * Sets the node color.
	 *
	 * @param color the node color
	 * @param pos the node position
	 */
	public void setNodeColor(final Color color, final int pos) {
		getNodeWithChecks(pos).setColor(color);
	}

	/**
	 * Sets the node color.
	 *
	 * @param color the node color
	 * @param pos the node position
	 */
	public void setNodeColor(final ColorRGB color, final int pos) {
		setNodeColor((color == null) ? null : new Color(color.getARGB()), pos);
	}

	/**
	 * Assigns the "value" property to this node.
	 *
	 * @param value the node value
	 * @param pos the node position
	 * @see PathProfiler#assignValues()
	 */
	public void setNodeValue(final double value, final int pos) {
		getNodeWithChecks(pos).v = value;
	}

	/**
	 * Assigns an annotation to an existing node.
	 *
	 * @param annotation the node annotation
	 * @param pos the node position
	 */
	public void setNodeAnnotation(final BrainAnnotation annotation, final int pos) {
		getNodeWithChecks(pos).setAnnotation(annotation);
	}

	/**
	 * Assigns a hemisphere to an existing node.
	 *
	 * @param hemisphereFlag the node hemisphere flag.
	 * @param pos the node position
	 */
	public void setNodeHemisphere(final char hemisphereFlag, final int pos) {
		getNodeWithChecks(pos).setHemisphere(hemisphereFlag);
	}

	/**
	 * Returns the "value" property of this node.
	 *
	 * @param pos the node position
	 * @return the value property of this node
	 * @see PathProfiler#assignValues()
	 * @see PointInImage#v
	 */
	public double getNodeValue(final int pos) {
		return getNodeWithChecks(pos).v;
	}

	/**
	 * Returns the neuropil annotation associated with this node.
	 *
	 * @param pos the node position
	 * @return the annotation of this node,
	 * @see SNTPoint#getAnnotation()
	 */
	public BrainAnnotation getNodeAnnotation(final int pos) {
		return getNodeWithChecks(pos).getAnnotation();
	}

	/**
	 * Returns the hemisphere flag associated with this node.
	 *
	 * @param pos the node position
	 * @return the flag associated with this node, either
	 *         {@link BrainAnnotation#LEFT_HEMISPHERE},
	 *         {@link BrainAnnotation#RIGHT_HEMISPHERE}, or
	 *         {@link BrainAnnotation#ANY_HEMISPHERE}
	 * @see SNTPoint#getAnnotation()
	 */
	public char getNodeHemisphereFlag(final int pos) {
		return getNodeWithChecks(pos).getHemisphere();
	}

	/**
	 * Sets the path "values", the array containing the numeric property assigned
	 * to path nodes, typically voxel intensities.
	 *
	 * @param values the new node values
	 * @see PointInImage#v
	 * @see PathProfiler#assignValues()
	 */
	public void setNodeValues(final double[] values) {
		if (values != null && values.length != size()) {
			throw new IllegalArgumentException("values array must have as many elements as nodes");
		}
		if (values == null) { // Clear all node values
			for (final PathNode node : nodes) node.v = Double.NaN;
		} else {
			for (int i = 0; i < values.length; i++) { // Set values on individual nodes
				getNodeWithoutChecks(i).v = values[i];
			}
		}
	}

	/**
	 * Assesses whether the nodes of this path have been assigned an array of
	 * values
	 *
	 * @return true, if successful
	 * @see PointInImage#v
	 * @see PathProfiler#assignValues()
	 */
	public boolean hasNodeValues() {
		// Direct iteration with early termination is faster than streams!?
		for (final PathNode node : nodes) {
			if (!Double.isNaN(node.v)) return true;
		}
		return false;
	}

	/**
	 * Assesses whether the nodes of this path have been assigned
	 * {@link BrainAnnotation}s.
	 *
	 * @return true, if successful
	 * @see SNTPoint#getAnnotation()
	 */
	public boolean hasNodeAnnotations() {
		// Direct iteration with early termination is faster than streams!?
		for (final PathNode node : nodes) {
			if (node.getAnnotation() != null) return true;
		}
		return false;
	}

	/**
	 * Assesses whether the nodes of this path have been assigned hemisphere flags.
	 *
	 * @return true, if successful
	 * @see SNTPoint#getHemisphere()
	 */
	public boolean hasNodeHemisphereFlags() {
		// Direct iteration with early termination is faster than streams!?
		for (final PathNode node : nodes) {
			if (node.getHemisphere() != BrainAnnotation.ANY_HEMISPHERE)
				return true;
		}
		return false;
	}

	/**
	 * Gets the color of this Path
	 *
	 * @return the color, or null if no color has been assigned to this Path
	 * @see #hasCustomColor
	 * @see #hasNodeColors()
	 */
	public Color getColor() {
		return color;
	}

	/**
	 * Gets the color of this Path
	 *
	 * @return the color, or null if no color has been assigned to this Path
	 * @see #hasCustomColor
	 * @see #hasNodeColors()
	 */
	public ColorRGB getColorRGB() {
		return (color == null) ? null : new ColorRGB(color.getRed(), color.getGreen(), color.getBlue());
	}

	/**
	 * Sets this path color.
	 *
	 * @param color the path color. Set it to null, to have SNT render this Path
	 *              using default settings. You may need to cast the null reference
	 *              to {@code Color} to avoid ambiguous method overload
	 */
	public void setColor(final Color color) {
		this.color = color;
		hasCustomColor = color != null;
		if (getFitted() != null) getFitted().setColor(color);
		if (hasNodeColors() && size() == 1) getNodeWithoutChecks(0).setColor(color);
	}

	/**
	 * Sets this path color.
	 *
	 * @param color the path color. Set it to null, to have SNT render this Path
	 *              using default settings. You may need to cast the null reference
	 *              to {@code ColorRGB} to avoid ambiguous method overload
	 */
	public void setColor(final ColorRGB color) {
		if (color == null)
			setColor((Color)null);
		else if (color instanceof ColorRGBA)
			setColor(new Color(color.getRed(), color.getGreen(), color.getBlue(), color.getAlpha()));
		else
			setColor(new Color(color.getARGB()));
	}

	/**
	 * Resets this path color, forcing SNT to render it under default settings.
	 */
	public void resetColor() {
		setColor((Color)null);
	}

	/**
	 * Assesses whether a custom color has been assigned to this Path.
	 *
	 * @return true, if successful
	 * @see #hasNodeColors()
	 */
	public boolean hasCustomColor() {
		return (hasCustomColor && color != null); // backwards compatibility: we cannot include hasNodeColors()
	}

	/**
	 * Gets the default SWC colors used by SNT.
	 *
	 * @param swcType the SEC type (e.g., {@link Path#SWC_AXON},
	 *          {@link Path#SWC_DENDRITE}, etc.)
	 * @return the SWC color
	 */
	public static Color getSWCcolor(final int swcType) {
        return switch (swcType) {
            case Path.SWC_SOMA -> Color.BLUE;
            case Path.SWC_AXON -> Color.RED;
            case Path.SWC_DENDRITE -> Color.GREEN;
            case Path.SWC_APICAL_DENDRITE -> Color.CYAN;
            case Path.SWC_CUSTOM -> Color.YELLOW;
            case Path.SWC_UNSPECIFIED -> Color.ORANGE;
            case Path.SWC_GLIA_PROCESS -> Color.PINK;
            case Path.SWC_CUSTOM2 -> Color.YELLOW.darker();
            default -> SNT.DEFAULT_DESELECTED_COLOR;
        };
	}

	/**
	 * Checks if this Path is a fitted version of another path.
	 *
	 * @return true, if it is a fitted version of another path
	 */
	public boolean isFittedVersionOfAnotherPath() {
		return fittedVersionOf != null;
	}

	protected void setFitted(final Path p) {
		if (getFitted() != null && p != null) {
			throw new IllegalArgumentException(
				"BUG: Trying to set a fitted path when there already is one...");
		}
		fitted = p;
		if (p == null) {
			setUseFitted(false);
		} else {
			p.fittedVersionOf = this;
			p.setIDs(-getID(), getTreeID());
//			System.out.println("Just set " + p.getName() + " is a fitted version of " + getName() );
		}
	}

	/**
	 * Sets whether the fitted flavor of this Path should be used.
	 *
	 * @param useFitted the new use fitted
	 * @throws IllegalArgumentException if this Path has not been fitted
	 */
	public void setUseFitted(final boolean useFitted)
		throws IllegalArgumentException
	{
		if (useFitted && getFitted() == null) throw new IllegalArgumentException(
			"setUseFitted(true) called, but 'fitted' member was null");
		this.useFitted = useFitted;
	}

	protected void discardFit() {
		fitted = null;
		useFitted = false;
	}

	/**
	 * Gets whether the fitted 'flavor' of this Path is in use.
	 *
	 * @return true, if the fitted version of this Path is in use, otherwise false
	 */
	public boolean getUseFitted() {
		return useFitted && fitted != null;
	}

	/**
	 * Gets the number of spines or varicosities associated with this path.
	 * @return the number of spines or varicosities. Returns 0 if no features have
	 *         been annotated.
	 * 
	 * @see #setSpineOrVaricosityCount(int) for setting the count
	 */
	public int getSpineOrVaricosityCount() {
		return spinesOrVaricosities;
	}

	/**
	 * Sets the number of spines or varicosities associated with this path.
	 * 
	 * @param newCount the new count of spines or varicosities. Should be non-negative.
	 * 
	 * @see #getSpineOrVaricosityCount() for getting the current count
	 */
	public void setSpineOrVaricosityCount(final int newCount) {
		this.spinesOrVaricosities = newCount;
		if (fitted != null) fitted.setSpineOrVaricosityCount(newCount);
	}

	/**
	 * Gets the hyperstack channel associated with this Path.
	 *
	 * @return the channel position of this path (1-based index). Note that if the
	 *         channel associated with a path is not known, it is assumed to be 1;
	 */
	public int getChannel() {
		return ctPosition[0];
	}

	/**
	 * Gets the hyperstack frame position associated with this Path.
	 *
	 * @return the frame position of this path (1-based index). Note that if the
	 *         frame associated with a path is not known, it is assumed to be 1;
	 */
	public int getFrame() {
		return ctPosition[1];
	}

	/**
	 * Sets the hyperstack position of this Path
	 *
	 * @param channel the channel (one-based index)
	 * @param frame the frame (one-based index)
	 */
	public void setCTposition(final int channel, final int frame) {
		ctPosition[0] = channel;
		ctPosition[1] = frame;
        if (fitted != null) {
            fitted.ctPosition[0] = channel;
            fitted.ctPosition[1] = frame;
        }
        if (fittedVersionOf != null) {
            fittedVersionOf.ctPosition[0] = channel;
            fittedVersionOf.ctPosition[1] = frame;
        }
	}

	/**
	 * Gets the fitted version ('flavor') of this Path.
	 *
	 * @return the fitted version, or null if this Path has not been fitted
	 */
	public Path getFitted() {
		return fitted;
	}

	public void setGuessedTangents(final int pointsEitherSide) {
		final double[] tangent = new double[3];
		for (int i = 0; i < nodes.size(); ++i) {
			getTangent(i, pointsEitherSide, tangent);
			nodes.get(i).setTangent(tangent.clone());
		}
	}

	public void getTangent(final int i, final int pointsEitherSide, final double[] result) {
		int min_index = i - pointsEitherSide;
		if (min_index < 0) min_index = 0;
		int max_index = i + pointsEitherSide;
		if (max_index >= size()) max_index = size() - 1;
		final PathNode minNode = nodes.get(min_index);
		final PathNode maxNode = nodes.get(max_index);
		result[0] = maxNode.x - minNode.x;
		result[1] = maxNode.y - minNode.y;
		result[2] = maxNode.z - minNode.z;
	}

	/**
	 * Gets the list of string representations of SWC types
	 *
	 * @return the list of SWC type names
	 */
	public static List<String> getSWCtypeNames() {
		return List.of(SWC_UNDEFINED_LABEL, SWC_SOMA_LABEL, SWC_AXON_LABEL, SWC_DENDRITE_LABEL,
				SWC_APICAL_DENDRITE_LABEL, SWC_CUSTOM_LABEL);
	}

	/**
	 * Gets the list of SWC types.
	 *
	 * @return the list of SWC type flags
	 */
	public static List<Integer> getSWCtypes() {
		return List.of(SWC_UNDEFINED, SWC_SOMA, SWC_AXON, SWC_DENDRITE, SWC_APICAL_DENDRITE, SWC_CUSTOM);
	}

	/**
	 * Gets the SWC type label associated with the specified type flag. SNT
	 * follows the specification detailed at <a href=
	 * "https://swc-specification.readthedocs.io/en/latest/">swc-specification.readthedocs.io</a>
	 *
	 * @param type the SWC type flag
	 * @param capitalized whether output String should be capitalized
	 * @return the respective label, or {@link Path#SWC_UNDEFINED_LABEL} if flag
	 *         was not recognized
	 */
	public static String getSWCtypeName(final int type,
		final boolean capitalized)
	{
		String typeName;
        switch (type) {
            case -1 -> {
                return (capitalized) ? "None" : "none";
            }
            case SWC_SOMA -> typeName = SWC_SOMA_LABEL;
            case SWC_AXON -> typeName = SWC_AXON_LABEL;
            case SWC_DENDRITE -> typeName = SWC_DENDRITE_LABEL;
            case SWC_APICAL_DENDRITE -> typeName = SWC_APICAL_DENDRITE_LABEL;
            case SWC_CUSTOM -> typeName = SWC_CUSTOM_LABEL;
            case SWC_UNSPECIFIED -> typeName = SWC_UNSPECIFIED_LABEL;
            case SWC_GLIA_PROCESS -> typeName = SWC_GLIA_PROCESS_LABEL;
            case SWC_CUSTOM2 -> typeName = SWC_CUSTOM2_LABEL;
            default -> typeName = SWC_UNDEFINED_LABEL;
        }
		return (capitalized) ? org.apache.commons.text.WordUtils.capitalize(typeName) : typeName;
	}

	/**
	 * Gets the path mean radius.
	 *
	 * @return the average radius of the path, or zero if path has no defined
	 *         thickness
	 * @see #hasRadii()
	 */
	public double getMeanRadius() {
		if (!hasRadii()) return 0;
		// Direct calculation without intermediate array creation
		double sum = 0;
		int count = 0;
		for (final PathNode node : nodes) {
			if (node.hasRadius()) {
				sum += node.getRadius();
				count++;
			}
		}
		return count > 0 ? sum / count : 0;
	}

	/**
	 * Gets the radius of the specified node.
	 *
	 * @param pos the node position
	 * @return the radius at the specified position, or zero if path has no
	 *         defined thickness
	 */
	public double getNodeRadius(final int pos) {
		return getNodeWithChecks(pos).getRadius();
	}

	/**
	 * Sets the radius of a specific node.
	 *
	 * @param radius the radius value
	 * @param pos the node position
	 */
	public void setNodeRadius(final double radius, final int pos) {
		getNodeWithChecks(pos).setRadius(radius);
	}

	/**
	 * Checks whether the nodes of this path have been assigned defined thickness.
	 *
	 * @return true, if the points defining with this path are associated with a
	 *         list of radii
	 */
	public boolean hasRadii() {
		// Direct iteration with early termination is faster than streams
		for (final PathNode node : nodes) {
			if (node.hasRadius())  return true;
		}
		return false;
	}

	/**
	 * Checks whether any nodes have custom colors assigned.
	 *
	 * @return true if any nodes have custom colors
	 * @see #getNodeColors()
	 */
	public boolean hasNodeColors() {
		// Direct iteration with early termination is faster than streams!?
		for (final PathNode node : nodes) {
			if (node.hasColor()) return true;
		}
		return false;
	}

	/**
	 * Gets the tangent vector of a specific node.
	 *
	 * @param pos the node position
	 * @return the tangent vector [x, y, z] or null if not defined
	 */
	public double[] getNodeTangent(final int pos) {
		return getNodeWithChecks(pos).getTangent();
	}

	/**
	 * Sets the tangent vector of a specific node.
	 *
	 * @param tangent the tangent vector [x, y, z]
	 * @param pos the node position
	 */
	public void setNodeTangent(final double[] tangent, final int pos) {
		getNodeWithChecks(pos).setTangent(tangent);
	}

	/**
	 * Sets the tangent vector of a specific node.
	 *
	 * @param tx the X component of the tangent
	 * @param ty the Y component of the tangent
	 * @param tz the Z component of the tangent
	 * @param pos the node position
	 */
	public void setNodeTangent(final double tx, final double ty, final double tz, final int pos) {
		getNodeWithChecks(pos).setTangent(tx, ty, tz);
	}

	/**
	 * Checks whether any nodes have tangent vectors defined.
	 *
	 * @return true if any nodes have tangents
	 */
	public boolean hasTangents() {
		// Direct iteration with early termination is faster than streams
		for (final PathNode node : nodes) {
			if (node.hasTangent()) return true;
		}
		return false;
	}

	/**
	 * Returns the first node of this path.
	 *
	 * @return the root node, or null if path is empty
	 */
	public PointInImage firstNode() {
		return size() > 0 ? getNodeWithoutChecks(0) : null;
	}

	/**
	 * Returns the last node of this path.
	 *
	 * @return the terminal node, or null if path is empty
	 */
	public PointInImage lastNode() {
		return size() > 0 ? getNodeWithoutChecks(size() - 1) : null;
	}

	/**
	 * Checks if this path contains the specified point within the given tolerance.
	 *
	 * @param point the point to check
	 * @param tolerance the distance tolerance
	 * @return true if the path contains the point within tolerance
	 */
	public boolean contains(final PointInImage point, final double tolerance) {
		for (int i = 0; i < size(); i++) {
			if (getNodeWithoutChecks(i).distanceTo(point) <= tolerance) {
				return true;
			}
		}
		return false;
	}

	/**
	 * Returns a new Path in which nodes have been interpolated to ensure that
	 * the distance between consecutive nodes does not exceed the specified
	 * threshold.
	 *
	 * @param maxDistance the maximum distance allowed between consecutive nodes
	 * @return the interpolated path
	 */
	public Path getInterpolatedPath(final double maxDistance) {
		if (maxDistance <= 0) {
			throw new IllegalArgumentException("maxDistance must be positive");
		}

		final Path result = createPath();
		if (size() == 0) return result;

		// Add first point
		result.addNode(getNodeWithoutChecks(0));

		for (int i = 1; i < size(); i++) {
			final PointInImage p1 = getNodeWithoutChecks(i - 1);
			final PointInImage p2 = getNodeWithoutChecks(i);

			final double distance = p1.distanceTo(p2);
			if (distance <= maxDistance) {
				// Distance is acceptable, add the point directly
				result.addNode(p2);
			} else {
				// Need to interpolate
				final int numSegments = (int) Math.ceil(distance / maxDistance);
				final double stepX = (p2.x - p1.x) / numSegments;
				final double stepY = (p2.y - p1.y) / numSegments;
				final double stepZ = (p2.z - p1.z) / numSegments;

				// Add interpolated points
				for (int j = 1; j <= numSegments; j++) {
					final double newX = p1.x + stepX * j;
					final double newY = p1.y + stepY * j;
					final double newZ = p1.z + stepZ * j;
					result.addNode(new PointInImage(newX, newY, newZ));
				}
			}
		}

		return result;
	}

	/**
	 * Returns a new Path in which nodes have been resampled to ensure uniform
	 * spacing along the path.
	 *
	 * @param targetSpacing the desired spacing between consecutive nodes
	 * @return the resampled path
	 */
	public Path getResampledPath(final double targetSpacing) {
		if (targetSpacing <= 0) {
			throw new IllegalArgumentException("targetSpacing must be positive");
		}

		final Path result = createPath();
		if (size() <= 1) {
			// Copy all nodes if path is too short
			for (int i = 0; i < size(); i++) {
				result.addNode(getNodeWithoutChecks(i));
			}
			return result;
		}

		// Calculate cumulative distances
		final double[] cumDistances = new double[size()];
		cumDistances[0] = 0;
		for (int i = 1; i < size(); i++) {
			final PointInImage p1 = getNodeWithoutChecks(i - 1);
			final PointInImage p2 = getNodeWithoutChecks(i);
			cumDistances[i] = cumDistances[i - 1] + p1.distanceTo(p2);
		}

		final double totalLength = cumDistances[size() - 1];
		if (totalLength == 0) {
			// All points are the same, just copy them
			for (int i = 0; i < size(); i++) {
				result.addNode(getNodeWithoutChecks(i));
			}
			return result;
		}

		// Add first point
		result.addNode(getNodeWithoutChecks(0));

		// Add resampled points
		double currentDistance = targetSpacing;
		while (currentDistance < totalLength) {
			// Find the segment containing this distance
			int segmentIndex = 0;
			for (int i = 1; i < size(); i++) {
				if (cumDistances[i] >= currentDistance) {
					segmentIndex = i - 1;
					break;
				}
			}

			// Interpolate within the segment
			final PointInImage p1 = getNodeWithoutChecks(segmentIndex);
			final PointInImage p2 = getNodeWithoutChecks(segmentIndex + 1);
			final double segmentStart = cumDistances[segmentIndex];
			final double segmentEnd = cumDistances[segmentIndex + 1];
			final double segmentLength = segmentEnd - segmentStart;

			if (segmentLength > 0) {
				final double t = (currentDistance - segmentStart) / segmentLength;
				final double newX = p1.x + t * (p2.x - p1.x);
				final double newY = p1.y + t * (p2.y - p1.y);
				final double newZ = p1.z + t * (p2.z - p1.z);
				result.addNode(new PointInImage(newX, newY, newZ));
			}

			currentDistance += targetSpacing;
		}

		// Add final point if it's not too close to the last added point
		final PointInImage lastOriginal = getNodeWithoutChecks(size() - 1);
		if (result.size() == 0 || result.getNodeWithoutChecks(result.size() - 1).distanceTo(lastOriginal) > targetSpacing / 2) {
			result.addNode(lastOriginal);
		}

		return result;
	}

	/**
	 * Returns a smoothed version of this path using a simple moving average filter.
	 *
	 * @param windowSize the size of the smoothing window (must be odd and >= 3)
	 * @return the smoothed path
	 */
	public Path getSmoothedPath(final int windowSize) {
		if (windowSize < 3 || windowSize % 2 == 0) {
			throw new IllegalArgumentException("windowSize must be odd and >= 3");
		}

		final Path result = createPath();
		if (size() <= windowSize) {
			// Path too short for smoothing, return copy
			for (int i = 0; i < size(); i++) {
				result.addNode(getNodeWithoutChecks(i));
			}
			return result;
		}

		final int halfWindow = windowSize / 2;

		// Copy first few points unchanged
		for (int i = 0; i < halfWindow; i++) {
			result.addNode(getNodeWithoutChecks(i));
		}

		// Apply smoothing to middle points
		for (int i = halfWindow; i < size() - halfWindow; i++) {
			double sumX = 0, sumY = 0, sumZ = 0;
			for (int j = i - halfWindow; j <= i + halfWindow; j++) {
				final PointInImage p = getNodeWithoutChecks(j);
				sumX += p.x;
				sumY += p.y;
				sumZ += p.z;
			}
			final double avgX = sumX / windowSize;
			final double avgY = sumY / windowSize;
			final double avgZ = sumZ / windowSize;
			result.addNode(new PointInImage(avgX, avgY, avgZ));
		}

		// Copy last few points unchanged
		for (int i = size() - halfWindow; i < size(); i++) {
			result.addNode(getNodeWithoutChecks(i));
		}

		return result;
	}

    /**
     * @deprecated Use {@link #sanitizeRadii(DoublePredicate, boolean)} instead
     */
    @Deprecated
    public Map<Integer, Double> interpolateMissingRadii(final DoublePredicate predicate,
                                                        final boolean apply) {
        return sanitizeRadii(predicate, apply);
    }

    /**
     * Sanitizes node radii by interpolating invalid values (zero, NaN, or negative).
     * <p>
     * Convenience method equivalent to calling
     * {@code sanitizeRadii((r) -> r <= 0 || Double.isNaN(r), apply)}.
     * </p>
     *
     * @param apply if {@code true}, corrected values are immediately applied
     * @return map of corrected {@code (node index, radius)} pairs, or {@code null}
     * if path has no radii or fewer than 2 nodes
     * @see #sanitizeRadii(DoublePredicate, boolean)
     */
    public Map<Integer, Double> sanitizeRadii(final boolean apply) {
        return sanitizeRadii(r -> r <= 0 || Double.isNaN(r), apply);
    }

    /**
     * Sanitizes node radii by interpolating values that match the given predicate.
     * <p>
     * Collects nodes with invalid radii (as defined by predicate) and assigns new
     * values using linear interpolation based on remaining valid nodes. Interpolated
     * values are clamped to the observed range of valid radii to prevent extrapolation
     * artifacts.
     * </p>
     * <p>
     * If fewer than two valid radii exist, interpolation is not possible and an
     * empty map is returned. For automatic fallback handling, use
     * {@link #sanitizeRadii(double, double, boolean)} instead.
     * </p>
     *
     * @param predicate function returning {@code true} for invalid radii
     * @param apply     if {@code true}, corrected values are immediately applied
     * @return map of corrected {@code (node index, radius)} pairs, or {@code null}
     * if path has no radii or fewer than 2 nodes. Empty map if no invalid
     * radii found or interpolation not possible.
     * @see #sanitizeRadii(double, double, boolean)
     */
    public Map<Integer, Double> sanitizeRadii(final DoublePredicate predicate, final boolean apply) {
        if (!hasRadii() || size() < 2)
            return null;

        // Collect valid and invalid nodes
        final List<Integer> validIndices = new ArrayList<>();
        final List<Double> validRadii = new ArrayList<>();
        final List<Integer> invalidIndices = new ArrayList<>();

        for (int nodeIdx = 0; nodeIdx < size(); nodeIdx++) {
            final double r = getNodeWithoutChecks(nodeIdx).getRadius();
            if (predicate.test(r)) {
                invalidIndices.add(nodeIdx);
            } else {
                validIndices.add(nodeIdx);
                validRadii.add(r);
            }
        }

        // Nothing to fix
        if (invalidIndices.isEmpty()) {
            return Collections.emptyMap();
        }

        // Cannot interpolate without at least 2 valid points
        if (validIndices.size() < 2) {
            SNTUtils.warn("Cannot interpolate radii: fewer than 2 valid radii found. " +
                    "Consider using sanitizeRadii(minRadius, maxRadius, apply) for automatic fallback.");
            return Collections.emptyMap();
        }

        // Compute bounds from valid radii for clamping
        final double minValid = validRadii.stream().mapToDouble(Double::doubleValue).min().orElse(0);
        final double maxValid = validRadii.stream().mapToDouble(Double::doubleValue).max().orElse(Double.MAX_VALUE);

        return interpolateAndApply(validIndices, validRadii, invalidIndices, minValid, maxValid, apply);
    }

    /**
     * Sanitizes node radii by enforcing bounds and interpolating outliers.
     * <p>
     * Identifies nodes with radii outside {@code [minRadius, maxRadius]} and assigns
     * new values using linear interpolation based on remaining valid nodes. Interpolated
     * values are clamped to the specified range to prevent extrapolation artifacts.
     * </p>
     * If fewer than two nodes have valid radii, a context-appropriate fallback is used:
     * <ul>
     *   <li>All radii above range → {@code maxRadius}</li>
     *   <li>All radii below range → {@code minRadius}</li>
     *   <li>Mixed or NaN values → midpoint of range</li>
     * </ul>
     *
     * @param minRadius minimum valid radius (inclusive)
     * @param maxRadius maximum valid radius (inclusive)
     * @param apply     if {@code true}, corrected values are immediately applied
     * @return map of corrected {@code (node index, radius)} pairs, or {@code null}
     * if path has no radii or fewer than 2 nodes. Empty map if all radii
     * are already within range.
     * @throws IllegalArgumentException if bounds are negative or {@code minRadius > maxRadius}
     * @see #sanitizeRadii(DoublePredicate, boolean)
     */
    public Map<Integer, Double> sanitizeRadii(final double minRadius,
                                              final double maxRadius,
                                              final boolean apply) {
        // Validate arguments
        if (minRadius < 0 || maxRadius < 0) {
            throw new IllegalArgumentException("Radius bounds must be non-negative");
        }
        if (minRadius > maxRadius) {
            throw new IllegalArgumentException(
                    "minRadius (" + minRadius + ") must be <= maxRadius (" + maxRadius + ")");
        }

        if (!hasRadii() || size() < 2)
            return null;

        // Collect valid and invalid nodes based on range
        final List<Integer> validIndices = new ArrayList<>();
        final List<Double> validRadii = new ArrayList<>();
        final List<Integer> invalidIndices = new ArrayList<>();

        for (int nodeIdx = 0; nodeIdx < size(); nodeIdx++) {
            final double r = getNodeWithoutChecks(nodeIdx).getRadius();
            if (r < minRadius || r > maxRadius || Double.isNaN(r)) {
                invalidIndices.add(nodeIdx);
            } else {
                validIndices.add(nodeIdx);
                validRadii.add(r);
            }
        }

        // Nothing to fix
        if (invalidIndices.isEmpty()) {
            return Collections.emptyMap();
        }

        // Fallback when interpolation not possible
        if (validIndices.size() < 2) {
            return applyFallback(invalidIndices, minRadius, maxRadius, apply);
        }

        return interpolateAndApply(validIndices, validRadii, invalidIndices, minRadius, maxRadius, apply);
    }

    /**
     * Determines context-appropriate fallback and applies to all invalid nodes.
     */
    private Map<Integer, Double> applyFallback(final List<Integer> invalidIndices,
                                               final double minRadius,
                                               final double maxRadius,
                                               final boolean apply) {
        // Analyze invalid radii to determine best fallback
        boolean allAboveMax = true;
        boolean allBelowMin = true;

        for (final int nodeIdx : invalidIndices) {
            final double r = getNodeWithoutChecks(nodeIdx).getRadius();
            if (!Double.isNaN(r)) {
                if (r <= maxRadius) allAboveMax = false;
                if (r >= minRadius) allBelowMin = false;
            }
        }

        final double fallback;
        if (allAboveMax) {
            fallback = maxRadius;
        } else if (allBelowMin) {
            fallback = minRadius;
        } else {
            fallback = (minRadius + maxRadius) / 2.0;
        }

        SNTUtils.log("Fewer than 2 valid radii; using fallback value: " + fallback);

        final Map<Integer, Double> result = new TreeMap<>();
        for (final int nodeIdx : invalidIndices) {
            result.put(nodeIdx, fallback);
            if (apply) {
                setRadius(fallback, nodeIdx);
            }
        }
        return result;
    }

    /**
     * Performs interpolation and applies clamped results.
     */
    private Map<Integer, Double> interpolateAndApply(final List<Integer> validIndices,
                                                     final List<Double> validRadii,
                                                     final List<Integer> invalidIndices,
                                                     final double clampMin,
                                                     final double clampMax,
                                                     final boolean apply) {
        // Build arrays for interpolation
        final double[] knownIndices = new double[validIndices.size()];
        final double[] knownRadii = new double[validRadii.size()];
        final double[] unknownIndices = new double[invalidIndices.size()];

        for (int i = 0; i < validIndices.size(); i++) {
            knownIndices[i] = validIndices.get(i);
            knownRadii[i] = validRadii.get(i);
        }
        for (int i = 0; i < invalidIndices.size(); i++) {
            unknownIndices[i] = invalidIndices.get(i);
        }

        // Interpolate
        final double[] interpolatedRadii = interpolate(knownIndices, knownRadii, unknownIndices);

        // Apply results with clamping
        final Map<Integer, Double> result = new TreeMap<>();
        for (int i = 0; i < unknownIndices.length; i++) {
            final int nodeIdx = (int) unknownIndices[i];
            double r = interpolatedRadii[i];

            // Skip invalid interpolation results
            if (Double.isNaN(r) || r < 0) {
                continue;
            }

            // Clamp to valid range
            r = Math.max(clampMin, Math.min(clampMax, r));

            result.put(nodeIdx, r);
            if (apply) {
                setRadius(r, nodeIdx);
            }
        }

        return result;
    }

    private double[] interpolate(final double[] x1, final double[] y1, final double[] x2) {
        // see https://stackoverflow.com/a/73716167
        final PolynomialSplineFunction function = new LinearInterpolator().interpolate(x1, y1);
        final PolynomialFunction[] splines = function.getPolynomials();
        final PolynomialFunction firstFunction = splines[0];
        final PolynomialFunction lastFunction = splines[splines.length - 1];
        final double[] knots = function.getKnots();
        final double firstKnot = knots[0];
        final double lastKnot = knots[knots.length - 1];
        return Arrays.stream(x2).map(d -> {
            if (d > lastKnot) {
                return lastFunction.value(d - knots[knots.length - 2]);
            } else if (d < firstKnot) {
                return firstFunction.value(d - knots[0]);
            }
            return function.value(d);
        }).toArray();
    }

	protected void setFittedCircles(final int nPoints, final double[] tangents_x,
		final double[] tangents_y, final double[] tangents_z,
		final double[] radii, final double[] optimized_x,
		final double[] optimized_y, final double[] optimized_z)
	{
		// Clear existing nodes and create new ones from the fitted data
		nodes.clear();
		
		for (int i = 0; i < nPoints; i++) {
			PathNode node = new PathNode(optimized_x[i], optimized_y[i], optimized_z[i]);
			node.onPath = this;
			
			// Set radius if available
			if (radii != null && i < radii.length) {
				node.setRadius(radii[i]);
			}
			
			// Set tangent if available
			if (tangents_x != null && tangents_y != null && tangents_z != null && 
				i < tangents_x.length && i < tangents_y.length && i < tangents_z.length) {
				node.setTangent(tangents_x[i], tangents_y[i], tangents_z[i]);
			}
			
			nodes.add(node);
		}
	}

	/**
	 * This toString() method shows details of the path which is actually being
	 * displayed, not necessarily this path object.
	 *
	 * @return the string
	 */
	@Override
	public String toString() {
		final StringBuilder sb = new StringBuilder(getName());
		if (getUseFitted() && sb.lastIndexOf("(Fitted)") < 0)
			sb.append(" (Fitted)");
		if (size() == 1)
			sb.append(" [Single Point]");
		if (swcType != SWC_UNDEFINED)
			sb.append(" [").append(getSWCtypeName(swcType, false)).append("]");
		return sb.toString();
	}

	/**
	 * Sets the SWC type.
	 *
	 * @param newSWCType the new SWC type
	 */
	public void setSWCType(final int newSWCType) {
		setSWCType(newSWCType, true);
	}

	protected void setSWCType(final int newSWCType,
		final boolean alsoSetInFittedVersion)
	{
		swcType = (newSWCType < 0 ) ? SWC_UNDEFINED : newSWCType;
		if (alsoSetInFittedVersion) {
			/*
			 * If we've been asked to also set the fitted version, this should only be
			 * called on the non-fitted version of the path, so raise an error if it's been
			 * called on the fitted version by mistake instead:
			 */
			if (isFittedVersionOfAnotherPath() && fittedVersionOf
				.getSWCType() != newSWCType) throw new IllegalArgumentException(
					"BUG: only call setSWCType on the unfitted path");
			if (getFitted() != null) getFitted().setSWCType(newSWCType);
		}
	}

	/**
	 * Gets the SWC type.
	 *
	 * @return the SWC type
	 */
	public int getSWCType() {
		return swcType;
	}

	/*
	 * @Override public String toString() { int n = size(); String result = ""; if(
	 * name != null ) result += "\"" + name + "\" "; result += n + " points"; if( n
	 * > 0 ) { result += " from " + x_positions[0] + ", " + y_positions[0] + ", " +
	 * z_positions[0]; result += " to " + x_positions[n-1] + ", " + y_positions[n-1]
	 * + ", " + z_positions[n-1]; } return result; }
	 */

	/**
	 * Gets the "bifurcation" (branching) order of this Path. If registered in the
	 * GUI, this would correspond to the level of this Path in
	 * {@link PathManagerUI}'s JTree: E.g., a Path connected to a primary Path
	 * (order 1) will be assigned order 2, etc. N.B.: Albeit related to reverse
	 * Horton-Strahler classification, Path ordering is formally distinct, as it
	 * classifies <i>Paths</i> instead of <i>branches</i>.
	 *
	 * @return the order of this path. A primary path is always of order 1.
	 * @see sc.fiji.snt.analysis.StrahlerAnalyzer
	 */
	public int getOrder() {
		return order;
	}

	public void setOrder(final int order) {
		this.order = order;
		if (getFitted() != null) getFitted().order = order;
		else if (fittedVersionOf != null) fittedVersionOf.order = order;
	}

	/*
	 * These are various fields that have the current 3D representations of this
	 * path. They should only be updated by synchronized methods, currently:
	 *
	 * updateContent3D addTo3DViewer removeFrom3DViewer
	 */
	int paths3DDisplay = 1;
	Content content3D;
	Content content3DExtra;
	ImagePlus content3DMultiColored;
	ImagePlus content3DExtraMultiColored;
	String nameWhenAddedToViewer;
	String nameWhenAddedToViewerExtra;

	synchronized void removeIncludingFittedFrom3DViewer(
		final Image3DUniverse univ)
	{
		removeFrom3DViewer(univ);
		if (useFitted) getFitted().removeFrom3DViewer(univ);
	}

	synchronized void updateContent3D(final Image3DUniverse univ,
		final boolean visible, final int paths3DDisplay, final Color3f color,
		final ImagePlus colorImage)
	{

//		SNT.log("In updateContent3D, colorImage is: " + colorImage);
//		SNT.log("In updateContent3D, color is: " + color);

		// So, go through each of the reasons why we might
		// have to remove (and possibly add back) the path:

		if (!visible) {
			/*
			 * It shouldn't be visible - if any of the contents are non-null, remove them:
			 */
			removeIncludingFittedFrom3DViewer(univ);
			return;
		}

		// Now we know it should be visible.
		Path pathToUse;

		if (useFitted) {
			/*
			 * If the non-fitted versions are currently being displayed, remove them:
			 */
			removeFrom3DViewer(univ);
			pathToUse = getFitted();
		}
		else {
			/*
			 * If the fitted version is currently being displayed, remove it:
			 */
			if (getFitted() != null) {
				getFitted().removeFrom3DViewer(univ);
			}
			pathToUse = this;
		}

//		if (SNT.isDebugMode()) {
//			SNT.log("pathToUse is: " + pathToUse);
//			SNT.log("  pathToUse.content3D is: " + pathToUse.content3D);
//			SNT.log("  pathToUse.content3DExtra is: " + pathToUse.content3DExtra);
//			SNT.log("  pathToUse.content3DMultiColored: " + pathToUse.content3DMultiColored);
//		}

		// Is the display (lines-and-discs or surfaces) right?
		if (pathToUse.paths3DDisplay != paths3DDisplay) {
			pathToUse.removeFrom3DViewer(univ);
			pathToUse.paths3DDisplay = paths3DDisplay;
			pathToUse.addTo3DViewer(univ, color, colorImage);
			return;
		}

		/* Were we previously using a colour image, but now not? */

		if (colorImage == null) {
			if ((paths3DDisplay == SNT.DISPLAY_PATHS_LINES_AND_DISCS &&
				pathToUse.content3DExtraMultiColored != null) ||
				(paths3DDisplay == SNT.DISPLAY_PATHS_SURFACE &&
					pathToUse.content3DMultiColored != null))
			{
				pathToUse.removeFrom3DViewer(univ);
				pathToUse.addTo3DViewer(univ, color, colorImage);
				return;
			}

			/*
			 * ... or, should we now use a colour image, where previously we were using a
			 * different colour image or no colour image?
			 */

		}
		else {
			if ((paths3DDisplay == SNT.DISPLAY_PATHS_LINES_AND_DISCS &&
				pathToUse.content3DExtraMultiColored != colorImage) ||
				(paths3DDisplay == SNT.DISPLAY_PATHS_SURFACE &&
					pathToUse.content3DMultiColored != colorImage))
			{
				pathToUse.removeFrom3DViewer(univ);
				pathToUse.addTo3DViewer(univ, color, colorImage);
				return;
			}
		}

		// Is the (flat) color wrong?
		if (pathToUse.realColor == null || !pathToUse.realColor.equals(color)) {

			/*
			 * If there's a representation of the path in the 3D viewer anyway, just set the
			 * color, don't recreate it, since the latter takes a long time:
			 */

			if (pathToUse.content3D != null || pathToUse.content3DExtra != null) {

				if (pathToUse.content3D != null) pathToUse.content3D.setColor(color);
				if (pathToUse.content3DExtra != null) pathToUse.content3DExtra.setColor(
					color);
				pathToUse.realColor = color;
				return;

			}
			// ... but if it wasn't in the 3D viewer, recreate it:
			pathToUse.removeFrom3DViewer(univ);
			pathToUse.paths3DDisplay = paths3DDisplay;
			pathToUse.addTo3DViewer(univ, color, colorImage);
			return;
		}

		if (pathToUse.nameWhenAddedToViewer == null || !univ.contains(
			pathToUse.nameWhenAddedToViewer))
		{
			pathToUse.paths3DDisplay = paths3DDisplay;
			pathToUse.addTo3DViewer(univ, color, colorImage);
		}
	}

	/**
	 * Removes this path from the specified 3D viewer universe.
	 * @param univ the Image3DUniverse from which to remove this path's 3D representation.
	 *             Must not be null.
	 *
	 * @deprecated This method is deprecated. Use newer 3D visualization methods instead.
	 */
	@Deprecated
	synchronized public void removeFrom3DViewer(final Image3DUniverse univ) {
		if (content3D != null) {
			univ.removeContent(nameWhenAddedToViewer);
			content3D = null;
		}
		if (content3DExtra != null) {
			univ.removeContent(nameWhenAddedToViewerExtra);
			content3DExtra = null;
		}
	}

	@Deprecated
	public java.util.List<Point3f> getPoint3fList() {
		final ArrayList<Point3f> linePoints = new ArrayList<>();
		for (final PathNode node : nodes) {
			linePoints.add(new Point3f((float) node.x, (float) node.y, (float) node.z));
		}
		return linePoints;
	}

	@Deprecated
	public Content addAsLinesTo3DViewer(final Image3DUniverse univ, final Color c,
		final ImagePlus colorImage)
	{
		return addAsLinesTo3DViewer(univ, Utils.toColor3f(c), colorImage);
	}

	@Deprecated
	protected Content addAsLinesTo3DViewer(final Image3DUniverse univ,
		final Color3f c, final ImagePlus colorImage)
	{
		final String safeName = univ.getSafeContentName(getName() + " as lines");
		return univ.addLineMesh(getPoint3fList(), c, safeName, true);
	}

	@Deprecated
	public Content addDiscsTo3DViewer(final Image3DUniverse univ, final Color c,
		final ImagePlus colorImage)
	{
		return addDiscsTo3DViewer(univ, Utils.toColor3f(c), colorImage);
	}

	@Deprecated
	public Content addDiscsTo3DViewer(final Image3DUniverse univ, final Color3f c,
		final ImagePlus colorImage)
	{
		if (!hasRadii()) return null;

		// Create coordinate arrays for color calculation - direct array creation is faster
		final int nodeCount = nodes.size();
		final double[] xCoords = new double[nodeCount];
		final double[] yCoords = new double[nodeCount];
		final double[] zCoords = new double[nodeCount];
		for (int i = 0; i < nodeCount; i++) {
			final PathNode node = nodes.get(i);
			xCoords[i] = node.x;
			yCoords[i] = node.y;
			zCoords[i] = node.z;
		}
		
		final Color3f[] originalColors = Pipe.getPointColors(xCoords, yCoords, zCoords, c, colorImage);

		final List<Color3f> meshColors = new ArrayList<>();

		final int edges = 8;
		final List<Point3f> allTriangles = new ArrayList<>(edges * size());
		for (int i = 0; i < size(); ++i) {
			final PathNode node = getNodeWithoutChecks(i);
			final double[] tangent = node.getTangent();
			
			// Use tangent if available, otherwise calculate from neighbors
			double tx = 0, ty = 0, tz = 1; // default tangent
			if (tangent != null && tangent.length >= 3) {
				tx = tangent[0];
				ty = tangent[1];
				tz = tangent[2];
			}
			
			final List<Point3f> discMesh = customnode.MeshMaker.createDisc(
				node.x, node.y, node.z, tx, ty, tz, node.getRadius(), 8);
			final int pointsInDiscMesh = discMesh.size();
			for (int j = 0; j < pointsInDiscMesh; ++j)
				meshColors.add(originalColors[i]);
			allTriangles.addAll(discMesh);
		}
		return univ.addTriangleMesh(allTriangles, meshColors, univ
			.getSafeContentName("Discs for path " + getName()));
	}

	/**
	 * Adds this path to the specified 3D viewer universe with the given color.
	 * 
	 * @param univ the Image3DUniverse to add this path to
	 * @param c the color to use for rendering this path
	 * @param colorImage the color image for texture mapping (can be null)
	 * 
	 * @deprecated This method is deprecated. Use newer 3D visualization methods instead.
	 */
	@Deprecated
	synchronized public void addTo3DViewer(final Image3DUniverse univ, final Color c, final ImagePlus colorImage) {
		addTo3DViewer(univ, Utils.toColor3f((c == null) ? SNT.DEFAULT_DESELECTED_COLOR : c), colorImage);
	}

	/**
	 * Adds this path to the specified 3D viewer universe with the given Color3f.
	 * 
	 * @param univ the Image3DUniverse to add this path to
	 * @param c the Color3f to use for rendering this path. Must not be null.
	 * @param colorImage the color image for texture mapping (can be null)
	 *
	 * @deprecated This method is deprecated. Use newer 3D visualization methods instead.
	 */
	@Deprecated
	synchronized public void addTo3DViewer(final Image3DUniverse univ, final Color3f c, final ImagePlus colorImage) {
		if (c == null) throw new IllegalArgumentException(
			"In addTo3DViewer, Color3f can no longer be null");

		realColor = c;

		if (size() <= 1) {
			content3D = null;
			content3DExtra = null;
			return;
		}

		if (paths3DDisplay == SNT.DISPLAY_PATHS_LINES ||
			paths3DDisplay == SNT.DISPLAY_PATHS_LINES_AND_DISCS)
		{
			content3D = addAsLinesTo3DViewer(univ, realColor, colorImage);
			content3D.setLocked(true);
			nameWhenAddedToViewer = content3D.getName();
			if (paths3DDisplay == SNT.DISPLAY_PATHS_LINES_AND_DISCS) {
				content3DExtra = addDiscsTo3DViewer(univ, realColor, colorImage);
				content3DExtraMultiColored = colorImage;
				if (content3DExtra == null) {
					nameWhenAddedToViewerExtra = null;
				}
				else {
					content3DExtra.setLocked(true);
					nameWhenAddedToViewerExtra = content3DExtra.getName();
				}
			}
			// univ.resetView();
			return;
		}

		int pointsToUse = -1;

		double[] x_points_d = new double[size()];
		double[] y_points_d = new double[size()];
		double[] z_points_d = new double[size()];
		double[] radiuses_d = new double[size()];

		if (hasRadii()) {
			int added = 0;
			final int samplingInterval = calculateSamplingInterval();
			int lastIndexAdded = -samplingInterval;
			for (int i = 0; i < size(); ++i) {
				if ((size() <= samplingInterval) || (i - lastIndexAdded >= samplingInterval)) {
					final PathNode node = getNodeWithoutChecks(i);
					x_points_d[added] = node.x;
					y_points_d[added] = node.y;
					z_points_d[added] = node.z;
					radiuses_d[added] = node.getRadius();
					lastIndexAdded = i;
					++added;
				}
			}
			pointsToUse = added;
		}
		else {
			for (int i = 0; i < size(); ++i) {
				final PathNode node = getNodeWithoutChecks(i);
				x_points_d[i] = node.x;
				y_points_d[i] = node.y;
				z_points_d[i] = node.z;
				radiuses_d[i] = getMinimumSeparation() * 2;
			}
			pointsToUse = size();
		}

		if (pointsToUse == 2) {
			// If there are only two points, then makeTube
			// fails, so interpolate:
			final double[] x_points_d_new = new double[3];
			final double[] y_points_d_new = new double[3];
			final double[] z_points_d_new = new double[3];
			final double[] radiuses_d_new = new double[3];

			x_points_d_new[0] = x_points_d[0];
			y_points_d_new[0] = y_points_d[0];
			z_points_d_new[0] = z_points_d[0];
			radiuses_d_new[0] = radiuses_d[0];

			x_points_d_new[1] = (x_points_d[0] + x_points_d[1]) / 2;
			y_points_d_new[1] = (y_points_d[0] + y_points_d[1]) / 2;
			z_points_d_new[1] = (z_points_d[0] + z_points_d[1]) / 2;
			radiuses_d_new[1] = (radiuses_d[0] + radiuses_d[1]) / 2;

			x_points_d_new[2] = x_points_d[1];
			y_points_d_new[2] = y_points_d[1];
			z_points_d_new[2] = z_points_d[1];
			radiuses_d_new[2] = radiuses_d[1];

			x_points_d = x_points_d_new;
			y_points_d = y_points_d_new;
			z_points_d = z_points_d_new;
			radiuses_d = radiuses_d_new;

			pointsToUse = 3;
		}

		final double[] x_points_d_trimmed = new double[pointsToUse];
		final double[] y_points_d_trimmed = new double[pointsToUse];
		final double[] z_points_d_trimmed = new double[pointsToUse];
		final double[] radiuses_d_trimmed = new double[pointsToUse];

		System.arraycopy(x_points_d, 0, x_points_d_trimmed, 0, pointsToUse);
		System.arraycopy(y_points_d, 0, y_points_d_trimmed, 0, pointsToUse);
		System.arraycopy(z_points_d, 0, z_points_d_trimmed, 0, pointsToUse);
		System.arraycopy(radiuses_d, 0, radiuses_d_trimmed, 0, pointsToUse);

		/*
		 * Work out whether to resample or not. I've found that the resampling is only
		 * really required in cases where the points are at adjacent voxels. So, work
		 * out the mean distance between all the points but in image co-ordinates - if
		 * there are points only at adjacent voxels this will be between 1 and sqrt(3)
		 * ~= 1.73. However, after the "fitting" process here, we might remove many of
		 * these points, so I'll say that we won't resample if the mean is rather higher
		 * - above 3. Hopefully this is a good compromise...
		 */

		double total_length_in_image_space = 0;
		for (int i = 1; i < pointsToUse; ++i) {
			final double x_diff = (x_points_d_trimmed[i] - x_points_d_trimmed[i -
				1]) / x_spacing;
			final double y_diff = (y_points_d_trimmed[i] - y_points_d_trimmed[i -
				1]) / y_spacing;
			final double z_diff = (z_points_d_trimmed[i] - z_points_d_trimmed[i -
				1]) / z_spacing;
			total_length_in_image_space += Math.sqrt(x_diff * x_diff + y_diff *
				y_diff + z_diff * z_diff);
		}
		final double mean_inter_point_distance_in_image_space =
			total_length_in_image_space / (pointsToUse - 1);
//		SNT.log("For path " + this + ", got mean_inter_point_distance_in_image_space: "
//				+ mean_inter_point_distance_in_image_space);
		final boolean resample = mean_inter_point_distance_in_image_space < 3;

//		SNT.log("... so" + (resample ? "" : " not") + " resampling");

		final ArrayList<Color3f> tubeColors = new ArrayList<>();

		final double[][][] allPoints = Pipe.makeTube(x_points_d_trimmed,
			y_points_d_trimmed, z_points_d_trimmed, radiuses_d_trimmed, resample ? 2
				: 1, // resample - 1 means just
			// "use mean distance
			// between points", 3 is
			// three times that,
			// etc.
			12, // "parallels" (12 means cross-sections are dodecagons)
			resample, // do_resample
			realColor, colorImage, tubeColors);

		if (allPoints == null) {
			content3D = null;
			content3DExtra = null;
			return;
		}

		// Make tube adds an extra point at the beginning and end:

		final List<Color3f> vertexColorList = new ArrayList<>();
		final List<Point3f> triangles = Pipe.generateTriangles(allPoints, 1, // scale
			tubeColors, vertexColorList);

		nameWhenAddedToViewer = univ.getSafeContentName(getName());
		// univ.resetView();
		content3D = univ.addTriangleMesh(triangles, vertexColorList,
			nameWhenAddedToViewer);
		content3D.setLocked(true);
		content3DMultiColored = colorImage;

		content3DExtra = null;
		nameWhenAddedToViewerExtra = null;

		// univ.resetView();
	}

	public void setSelected(final boolean newSelectedStatus) {
		selected = newSelectedStatus;
	}

	public boolean isSelected() {
		return selected;
	}

	/**
	 * Checks if this version of the path is currently being used for display/analysis.
	 * 
	 * @return true if this version is active, false otherwise
	 */
	public boolean isActiveVersion() {
		if (fittedVersionOf != null) return fittedVersionOf.useFitted;
		return !useFitted;
	}
	
	/**
	 * @deprecated Use {@link #isActiveVersion()} instead
	 */
	@Deprecated
	public boolean versionInUse() {
		return isActiveVersion();
	}

	/**
	 * Returns an estimated volume of this path.
	 * <p>
	 * The most accurate volume of each path segment would be the volume of a
	 * convex hull of two arbitrarily oriented and sized circles in space. This is
	 * tough to work out analytically, and this precision isn't really warranted
	 * given the errors introduced in the fitting process, the tracing in the
	 * first place, etc. So, this method produces an approximate volume assuming
	 * that the volume of each of these parts is that of a truncated cone
	 * (Frustum) , with circles of the same size (i.e., as if the circles had
	 * simply been reoriented to be parallel and have a common normal vector)
	 * </p>
	 * <p>
	 * For more accurate measurements of the volumes of a neuron, you should use
	 * the filling interface.
	 * </p>
	 *
	 * @return the approximate fitted volume (in physical units), or NaN if this
	 *         Path has no radii
	 * @see #hasRadii()
	 */
	public double getApproximatedVolume() {
		if (!hasRadii()) {
			return Double.NaN;
		}

		if (nodes.size() < 2) return 0;
		
		double totalVolume = 0;
		// Cache node references to avoid repeated method calls
		PathNode prevNode = nodes.get(0);
		for (int i = 1; i < nodes.size(); ++i) {
			final PathNode currentNode = nodes.get(i);
			
			final double xdiff = currentNode.x - prevNode.x;
			final double ydiff = currentNode.y - prevNode.y;
			final double zdiff = currentNode.z - prevNode.z;
			final double h = Math.sqrt(xdiff * xdiff + ydiff * ydiff + zdiff * zdiff);
			
			final double r1 = prevNode.getRadius();
			final double r2 = currentNode.getRadius();
			
			// See http://en.wikipedia.org/wiki/Frustum
			final double partVolume = (Math.PI * h * (r1 * r1 + r2 * r2 + r1 * r2)) /
				3.0;
			totalVolume += partVolume;
			prevNode = currentNode;
		}

		return totalVolume;
	}

	/**
	 * Returns an estimated surface area of this path, treating each inter-node
	 * segment as a frustum.
	 *
	 * @return the approximate surface area (in physical units), or 0 if this Path
	 *         has no radii
	 * @see #hasRadii()
	 * @see #getApproximatedVolume()
	 */
	public double getApproximatedSurface() {
		if (!hasRadii()) {
			return Double.NaN;
		}

		if (nodes.size() < 2) return 0;
		
		double totalSurface = 0;
		// Cache node references to avoid repeated method calls
		PathNode prevNode = nodes.get(0);
		for (int i = 1; i < nodes.size(); ++i) {
			final PathNode currentNode = nodes.get(i);
			
			final double xdiff = currentNode.x - prevNode.x;
			final double ydiff = currentNode.y - prevNode.y;
			final double zdiff = currentNode.z - prevNode.z;
			final double h = Math.sqrt(xdiff * xdiff + ydiff * ydiff + zdiff * zdiff);
			
			final double r1 = prevNode.getRadius();
			final double r2 = currentNode.getRadius();
			
			// Lateral surface area of conical frustum https://en.wikipedia.org/wiki/Frustum
			final double partSurface = Math.PI * (r1 + r2) * Math.sqrt((r1 - r2) * (r1 - r2) + (h * h));
			totalSurface += partSurface;
			prevNode = currentNode;
		}

		return totalSurface;
	}
	/*
	 * This doesn't deal with the startJoins, endJoins or fitted fields, since they
	 * involve other paths which were probably also transformed by the caller.
	 */

	public Path transform(final PathTransformer transformation,
		final ImagePlus template, final ImagePlus model)
	{

		double templatePixelWidth = 1;
		double templatePixelHeight = 1;
		double templatePixelDepth = 1;
		String templateUnits = "pixels";

		final Calibration templateCalibration = template.getCalibration();
		if (templateCalibration != null) {
			templatePixelWidth = templateCalibration.pixelWidth;
			templatePixelHeight = templateCalibration.pixelHeight;
			templatePixelDepth = templateCalibration.pixelDepth;
			templateUnits = SNTUtils.getSanitizedUnit(templateCalibration.getUnit());
		}

		final Path result = new Path(templatePixelWidth, templatePixelHeight,
			templatePixelDepth, templateUnits);
		final double[] transformed = new double[3];

		// Actually, just say you'll have to refit all the
		// previously fitted paths...

		for (int i = 0; i < size(); ++i) {
			final PathNode node = getNodeWithoutChecks(i);
			final double original_x = node.x;
			final double original_y = node.y;
			final double original_z = node.z;
			transformation.transformPoint(original_x, original_y, original_z,
				transformed);
			final double new_x = transformed[0];
			final double new_y = transformed[1];
			final double new_z = transformed[2];
			if (Double.isNaN(new_x) || Double.isNaN(new_y) || Double.isNaN(new_z))
				continue;
			result.addPointDouble(new_x, new_y, new_z);
		}

		result.id = id;
		result.selected = selected;
		result.name = name;

		result.x_spacing = x_spacing;
		result.y_spacing = y_spacing;
		result.z_spacing = z_spacing;
		result.spacing_units = spacing_units;

		result.swcType = swcType;

		return result;
	}

	/**
	 * Returns the nodes which are indicated to be a join (junction/branch point),
	 * either in this Path object, or any other that starts or ends on it.
	 *
	 * @return the list of nodes as {@link PointInImage} objects
	 * @see #findJunctionIndices()
	 * @see #getBranchPoints()
	 */
	public List<PointInImage> findJunctions() {
		final ArrayList<PointInImage> result = new ArrayList<>();
		final PointInImage myBranchPoint = getBranchPoint();
		if (myBranchPoint != null) {
			result.add(myBranchPoint);
		}
		for (final Path other : connectedPaths) {
			if (other.parentPath == this) {
				final PointInImage otherBranchPoint = other.getBranchPoint();
				if (otherBranchPoint != null) {
					result.add(otherBranchPoint);
				}
			}
		}
		return result;
	}

	/**
	 * This is a version of {@link #findJunctions()} ensuring that a junction node
	 * is only retrieved once even if multiple child paths are associated with it.
	 *
	 * @see #findJunctionIndices()
	 * @return the junction nodes
	 */
	public Set<PointInImage> getBranchPoints() {
		final Set<PathAndFillManager.PositionPoint> uniqueJunctions = new HashSet<>();
		final Set<PointInImage> result = new HashSet<>();
		for (final PointInImage junction : findJunctions()) {
			final PathAndFillManager.PositionPoint junctionPoint = new PathAndFillManager.PositionPoint(junction);
			if (uniqueJunctions.add(junctionPoint)) {
				result.add(junction); // This is a new unique junction based on location
			}
		}
		return result;
	}

	/** @deprecated use {@link #getBranchPoints()} instead */
	@Deprecated
	public Set<PointInImage> getJunctionNodes() {
		return getBranchPoints();
	}

	/**
	 * Returns the indices of nodes which are indicated to be a join, either in this
	 * Path object, or any other that starts or ends on it.
	 *
	 * @return the indices of junction nodes, naturally sorted
	 * @see #findJunctions()
	 * @see #getBranchPoints()
	 */
	public TreeSet<Integer> findJunctionIndices() {
		final TreeSet<Integer> result = new TreeSet<>();
		for (final PointInImage point : findJunctions()) {
			result.add(indexNearestTo(point.x, point.y, point.z));
		}
		return result;
	}

	/**
	 * Downsamples this path (in-place) by reducing the number of nodes while preserving its overall shape.
	 * <p>
	 * This method reduces the density of nodes in the path by removing redundant points expected
	 * to not significantly contribute to the path's shape. The downsampling is performed using the
	 * Douglas-Peucker algorithm, which preserves fidelity within the specified tolerance.
	 * </p>
	 * <p>
	 * The method operates on the segments flanked by anchor points (junctions, start, and end points).
	 * Each segment is downsampled independently, and the results are combined to form the final
	 * downsampled path. Node radii are averaged appropriately for retained points
	 * </p>
	 * <p>
	 * <b>Thread Safety:</b> This method is synchronized to prevent concurrent modification
	 * of the path structure during the downsampling operation.
	 * </p>
	 *
	 * @param internodeSpacing the target spacing between nodes after downsampling. This parameter
	 *                         controls the aggressiveness of the downsampling - smaller values
	 *                         preserve more detail, larger values result in more aggressive
	 *                         simplification. Must be greater than zero.
	 * @throws IllegalArgumentException if internodeSpacing is less than or equal to zero
	 * @see #upsample(double)
	 * @see PathDownsampler#downsample(ArrayList, double)
	 * @see #findJunctionIndices()
	 */
    public synchronized void downsample(final double internodeSpacing) {
        if (size() < 3) return; // Need at least 3 points to downsample

        // Find fixed points (junctions) that must be preserved. Set is already sorted
        // Convert to sorted array for easier processing
        final TreeSet<Integer> fixedIndicesSet = findJunctionIndices();
        fixedIndicesSet.add(0); // Always preserve start
        fixedIndicesSet.add(size() - 1); // Always preserve end
        // Direct array creation is faster than stream operations
        final int[] fixedIndices = new int[fixedIndicesSet.size()];
        int idx = 0;
        for (Integer index : fixedIndicesSet) {
            fixedIndices[idx++] = index;
        }

        // Process all segments and collect results
        final List<SegmentResult> segmentResults = new ArrayList<>();
        int totalPointsToKeep = 0;

        for (int segIdx = 0; segIdx < fixedIndices.length - 1; segIdx++) {
            final int segStart = fixedIndices[segIdx];
            final int segEnd = fixedIndices[segIdx + 1];

            // Skip empty or single-point segments
            if (segEnd <= segStart) continue;

            // Create input for downsampling
            final ArrayList<SimplePoint> segmentPoints = new ArrayList<>(segEnd - segStart + 1);
            for (int i = segStart; i <= segEnd; i++) {
                final PathNode node = getNodeWithoutChecks(i);
                segmentPoints.add(new SimplePoint(node.x, node.y, node.z, i));
            }

            // Downsample this segment
            ArrayList<SimplePoint> downsampled = PathDownsampler.downsample(segmentPoints, internodeSpacing);

            if (downsampled.isEmpty()) {
                // Fallback to original segment if downsampling fails
                downsampled = segmentPoints;
            }

            final SegmentResult result = new SegmentResult(segStart, segEnd, downsampled);
            segmentResults.add(result);
            totalPointsToKeep += downsampled.size();

            // Avoid double-counting shared endpoints between segments
            if (segIdx > 0) totalPointsToKeep--;
        }

        if (totalPointsToKeep >= size() || segmentResults.isEmpty()) return; // No improvement from downsampling

        // Create new nodes list
        final List<PathNode> newNodes = new ArrayList<>(totalPointsToKeep);

        // Build the final downsampled path
        for (int segIdx = 0; segIdx < segmentResults.size(); segIdx++) {
            final SegmentResult result = segmentResults.get(segIdx);
            final int startIdx = (segIdx == 0) ? 0 : 1; // Skip first point of non-first segments to avoid duplication

            for (int i = startIdx; i < result.downsampledPoints.size(); i++) {
                final SimplePoint sp = result.downsampledPoints.get(i);

                // Create new node with downsampled coordinates
                final PathNode newNode = new PathNode(sp.x, sp.y, sp.z);
                newNode.onPath = this;

                // Copy properties from original node if it exists
                if (sp.originalIndex >= 0 && sp.originalIndex < nodes.size()) {
                    final PathNode originalNode = nodes.get(sp.originalIndex);
                    newNode.v = originalNode.v;
                    newNode.setAnnotation(originalNode.getAnnotation());
                    newNode.setHemisphere(originalNode.getHemisphere());

                    // For radius, use average if this is an interpolated point
                    if (originalNode.hasRadius()) {
                        newNode.setRadius(calculateAverageRadius(sp, result, i));
                    }
                }

                newNodes.add(newNode);
            }
        }

        // Save children's branch points BEFORE modifying nodes
        final Map<Path, PointInImage> savedBranchPoints = new HashMap<>();
        for (final Path child : children) {
            final PointInImage bp = child.getBranchPoint();
            if (bp != null) {
                savedBranchPoints.put(child, bp.clone());
            }
        }

        // Replace nodes with downsampled version
        nodes.clear();
        nodes.addAll(newNodes);
        setNodeColors(null);

        // Recalculate children's branch point indices using saved coordinates
        for (final Path child : children) {
            final PointInImage savedBP = savedBranchPoints.get(child);
            if (savedBP != null) {
                child.branchPointIndex = indexNearestTo(savedBP.x, savedBP.y, savedBP.z);
            }
        }
    }

	/*
	 * Helper record to store the results of downsampling a single path segment.
	 * <p>
	 * This record encapsulates the original segment boundaries and the resulting
	 * downsampled points for a segment between two fixed points (junctions).
	 * </p>
	 */
	private record SegmentResult(int originalStart, int originalEnd, List<SimplePoint> downsampledPoints) {
	}
	
	/*
	 * Calculates an appropriate radius value for a downsampled point.
	 * <p>
	 * For endpoint nodes, the original radius is preserved. For intermediate nodes,
	 * the radius is computed as the average of nearby original nodes within a small
	 * neighborhood window to maintain smooth radius transitions.
	 * </p>
	 * 
	 * @param sp the downsampled point for which to calculate the radius
	 * @param result the segment result containing context about the original segment
	 * @param pointIndex the index of the point within the downsampled segment
	 * @return the calculated radius value, or 0.0 if no radius data is available
	 */
	private double calculateAverageRadius(final SimplePoint sp, final SegmentResult result, final int pointIndex) {
		if (!hasRadii() || sp.originalIndex < 0 || sp.originalIndex >= size()) {
			return 0.0;
		}
		
		// For endpoints, use the original radius
		if (pointIndex == 0 || pointIndex == result.downsampledPoints.size() - 1) {
			return getNodeWithoutChecks(sp.originalIndex).getRadius();
		}
		
		// For intermediate points, average over a small neighborhood
		final int originalIdx = sp.originalIndex;
		final int windowSize = 2; // Average over ±2 points
		final int startIdx = Math.max(0, Math.max(result.originalStart, originalIdx - windowSize));
		final int endIdx = Math.min(size() - 1, Math.min(result.originalEnd, originalIdx + windowSize));
		
		// Ensure valid range
		if (startIdx > endIdx) {
			return getNodeWithoutChecks(originalIdx).getRadius();
		}

		double sum = 0.0;
		int count = 0;
		for (int i = startIdx; i <= endIdx; i++) {
			sum += getNodeWithoutChecks(i).getRadius();
			count++;
		}
		return count > 0 ? sum / count : getNodeWithoutChecks(originalIdx).getRadius();
	}

	/**
	 * Upsamples this path by spacing nodes at a specified distance.
	 * <p>
	 * The original nodes are preserved and additional interpolated nodes are placed at regular
	 * intervals between them. The spacing between adjacent nodes will thus be approximately
	 * equal to the specified distance. If the distance between two adjacent original nodes is
	 * less than the specified pacing, no additional nodes are added between them.
	 * </p>
	 * <p>
	 *     <b>NB:</b>Some assigned properties (node colors, node values, etc.) will be lost.
	 * </p>
	 * <p>
	 * <b>Thread Safety:</b> This method is synchronized to prevent concurrent modification
	 * of the path structure during the upsampling operation.
	 * </p>
	 *
	 * @param internodeSpacing the desired distance between adjacent nodes
	 * @throws IllegalArgumentException if spacing is less than or equal to zero
	 */
	public synchronized void upsample(final double internodeSpacing) {
		if (internodeSpacing <= 0) {
			throw new IllegalArgumentException("internode spacing must be greater than zero");
		}

		if (size() < 2) return;

		// Prefer squared distances for (marginal) performance gains
		final double spacingSquared = internodeSpacing * internodeSpacing;
		final int originalSize = size();

		// Pre-calculate segment data and total points needed
		final double[] segmentDx = new double[originalSize - 1];
		final double[] segmentDy = new double[originalSize - 1];
		final double[] segmentDz = new double[originalSize - 1];
		final int[] segmentPoints = new int[originalSize - 1];
		int totalPointsToAdd = 0;

		for (int i = 0; i < originalSize - 1; i++) {
			final PathNode p1 = getNodeWithoutChecks(i);
			final PathNode p2 = getNodeWithoutChecks(i + 1);

			final double dx = p2.x - p1.x;
			final double dy = p2.y - p1.y;
			final double dz = p2.z - p1.z;
			final double lengthSquared = dx * dx + dy * dy + dz * dz;

			segmentDx[i] = dx;
			segmentDy[i] = dy;
			segmentDz[i] = dz;

			if (lengthSquared > spacingSquared) {
				final int numPoints = (int) Math.floor(Math.sqrt(lengthSquared) / internodeSpacing);
				segmentPoints[i] = numPoints;
				totalPointsToAdd += numPoints;
			} else {
				segmentPoints[i] = 0;
			}
		}

		if (totalPointsToAdd == 0) return;

		// Create new nodes list
		final List<PathNode> newNodes = new ArrayList<>(originalSize + totalPointsToAdd);

		// Process each segment
		for (int i = 0; i < originalSize - 1; i++) {
			final PathNode p1 = getNodeWithoutChecks(i);
			final PathNode p2 = getNodeWithoutChecks(i + 1);

			// Add original point
			PathNode newNode = p1.clone();
			newNode.onPath = this;
			newNodes.add(newNode);

			// Add interpolated points
			final int numPoints = segmentPoints[i];
			if (numPoints > 0) {
				final double stepSize = 1.0 / (numPoints + 1);
				final double dx = segmentDx[i];
				final double dy = segmentDy[i];
				final double dz = segmentDz[i];
				
				// Pre-calculate interpolation deltas
				final double r1 = p1.getRadius();
				final double r2 = p2.getRadius();
				final double v1 = p1.v;
				final double v2 = p2.v;
				final double dr = r2 - r1;
				final double dv = v2 - v1;
				
				// Pre-calculate midpoint for annotation/hemisphere selection
				final int midPoint = (numPoints + 1) / 2;

				for (int j = 1; j <= numPoints; j++) {
					final double t = j * stepSize;
					
					// Create interpolated node
					PathNode interpolatedNode = new PathNode(
						p1.x + t * dx,
						p1.y + t * dy,
						p1.z + t * dz
					);
					interpolatedNode.onPath = this;
					
					// Interpolate radius
					if (p1.hasRadius() || p2.hasRadius()) {
						interpolatedNode.setRadius(r1 + t * dr);
					}
					
					// Interpolate value
					interpolatedNode.v = v1 + t * dv;
					
					// Use closest node for annotation/hemisphere
					final PathNode closestNode = (j <= midPoint) ? p1 : p2;
					interpolatedNode.setAnnotation(closestNode.getAnnotation());
					interpolatedNode.setHemisphere(closestNode.getHemisphere());
					
					newNodes.add(interpolatedNode);
				}
			}
		}

		// Add final node
		PathNode lastNode = getNodeWithoutChecks(originalSize - 1).clone();
		lastNode.onPath = this;
		newNodes.add(lastNode);

		// Replace nodes with upsampled version
		nodes.clear();
		nodes.addAll(newNodes);
		setNodeColors(null);
	}

	private void setRadiiInPlace(final double[] newRadii) {
		if (newRadii != null) {
			assert size() == newRadii.length;
			// Set radii on individual nodes
			for (int i = 0; i < newRadii.length; i++) {
				getNodeWithoutChecks(i).setRadius(newRadii[i]);
			}
			// Update tangents if we now have radii
			setGuessedTangents(2);
		} else {
			// Clear all radii
			for (PathNode node : nodes) {
				node.setRadius(0.0);
			}
		}
	}

	/**
	 * Assigns a fixed radius to all the nodes of this Path
	 *
	 * @param r the radius to be assigned. Setting it to 0 or Double.NaN removes
	 *          the radius attribute from the Path
	 */
	public void setRadius(final double r) {
		if (Double.isNaN(r) || r == 0d) {
			// Clear all radii
			for (PathNode node : nodes) {
				node.setRadius(0.0);
			}
		}
		else {
			// Set radius for all nodes
			for (PathNode node : nodes) {
				node.setRadius(r);
			}
			// Update tangents if this is the first time setting radii
			if (!hasTangents()) {
				setGuessedTangents(2);
			}
		}
	}

	/**
	 * Assigns a radius to a specific node of this Path.
	 * <p>
	 * Sets the radius for the node at the specified index. If this is the first
	 * radius being set for the path and no tangents exist, this method will also
	 * compute guessed tangents to support proper 3D rendering.
	 * </p>
	 * 
	 * @param r the radius to be assigned to the node. NaN values are ignored.
	 * @param index the 0-based index of the node to modify
	 * 
	 * @throws IndexOutOfBoundsException if index is out of range
	 * 
	 * @see #setRadius(double) for setting radius on all nodes
	 * @see #getNodeRadius(int) for getting the radius of a specific node
	 */
	public void setRadius(final double r, final int index) {
		if (Double.isNaN(r)) return;
		getNodeWithChecks(index).setRadius(r);
		// Update tangents if this is the first radius being set
		if (!hasTangents()) setGuessedTangents(2);
	}

	/**
	 * Assigns radii to this Path
	 *
	 * @param radii the radii array. Setting it null removes the radius attribute
	 *          from the Path
	 * @see #setRadius(double)
	 */
	public void setRadii(final double[] radii) {
		if (radii == null || radii.length == 0) {
			for (final PathNode node : nodes)
				node.radius = 0; // Clear all radii
		} else if (radii.length != size()) {
			throw new IllegalArgumentException("radii array must have as many elements as nodes");
		} else {
			for (int i = 0; i < radii.length; i++) // Set radii on individual nodes
				getNodeWithoutChecks(i).setRadius(radii[i]);
			// Update tangents if we now have radii
			if (!hasTangents()) setGuessedTangents(2);
		}
	}

	/**
	 * Fits radii to this Path by fitting circular cross-sections to the image
	 * signal. This is the scripting equivalent of the UI's "Fit Path" action.
	 * Node positions are not modified; only radii and tangent vectors are
	 * computed.
	 *
	 * @param imp the image containing the signal to fit against
	 * @return true if the fit succeeded
	 * @throws IllegalArgumentException if this path is already a fitted version
	 *                                  of another path, or if the image is null
	 * @see PathFitter
	 * @see #setRadii(double[])
	 */
	public boolean fitRadii(final ImagePlus imp) {
		if (imp == null)
			throw new IllegalArgumentException("Image cannot be null");
		final PathFitter fitter = new PathFitter(imp, this);
		fitter.setScope(PathFitter.RADII);
		fitter.setReplaceNodes(true); // safe: radii-only, no midpoint changes
		fitter.call();
		if (fitter.getSucceeded()) {
			fitter.applyFit();
			return true;
		}
		return false;
	}

	public boolean isConnectedTo(final Path other) {
		return (getParentPath() != null && getParentPath().equals(other))
			|| (connectedPaths != null && connectedPaths.contains(other));
	}

	protected void addChangeListener(PathChangeListener listener) {
		changeListeners.add(listener);
	}

	protected boolean removeChangeListener(PathChangeListener listener) {
		return changeListeners.remove(listener);
	}

}