/FileExtraction/Java_unseen/fiji_AnalyzeSkeleton/src/main/java/sc/fiji/analyzeSkeleton/AnalyzeSkeleton_.java
/*
 * #%L
 * AnalyzeSkeleton_ plugin for ImageJ.
 * %%
 * Copyright (C) 2008 - 2017 Ignacio Arganda-Carreras.
 * %%
 * 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.analyzeSkeleton;

import ij.*;
import ij.gui.DialogListener;
import ij.gui.GenericDialog;
import ij.gui.Roi;
import ij.measure.Calibration;
import ij.measure.ResultsTable;
import ij.plugin.filter.PlugInFilter;
import ij.plugin.frame.Recorder;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;

import java.awt.AWTEvent;
import java.awt.Checkbox;
import java.awt.Font;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
import java.util.ListIterator;

/**
 * AnalyzeSkeleton_ plugin for ImageJ and Fiji.
 *
 * 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 (http://www.gnu.org/licenses/gpl.txt )
 *
 * 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 * 
 */

/**
 * Main class of the ImageJ/Fiji plugin for skeleton analysis.
 * This class is a plugin for the ImageJ and Fiji interfaces for analyzing
 * 2D/3D skeleton images.
 * <p>
 * For more detailed information, visit the AnalyzeSkeleton home page:
 * <A target="_blank" href="http://fiji.sc/AnalyzeSkeleton">http://fiji.sc/AnalyzeSkeleton</A>
 *
 *
 * @author Ignacio Arganda-Carreras
 *
 */
public class AnalyzeSkeleton_ implements PlugInFilter, DialogListener
{
	private static final String PRUNE_MODE_INDEX_KEY = "sc.fiji.analyzeSkeleton.pruneModeIndex";
	private static final String PRUNE_ENDS_KEY = "sc.fiji.analyzeSkeleton.pruneEnds";
	private static final String CALCULATE_PATH_KEY = "sc.fiji.analyzeSkeleton.shortestPath";
	private static final String VERBOSE_KEY = "sc.fiji.analyzeSkeleton.showDetailedInfo";
	private static final String DISPLAY_SKELETONS_KEY = "sc.fiji.analyzeSkeleton.displayLabeledSkeletons";
	private static final int DEFAULT_PRUNE_MODE_INDEX = AnalyzeSkeleton_.NONE;
	private static final boolean DEFAULT_PRUNE_ENDS = false;
	private static final boolean DEFAULT_CALCULATE_SHORTEST_PATH = false;
	private static final boolean DEFAULT_VERBOSE = false;
	private static final boolean DEFAULT_PROTECT_ROI = false;
	private static final boolean DEFAULT_DISPLAY_SKELETONS = false;
	private static final String HELP_URL = "http://fiji.sc/wiki/index.php/AnalyzeSkeleton";

	/** end point flag */
	public static byte END_POINT = 30;
	/** junction flag */
	public static byte JUNCTION = 70;
	/** slab flag */
	public static byte SLAB = 127;
	/** shortest path flag */
	public static byte SHORTEST_PATH = 96;

	private GenericDialog settingsDialog;

	/** working image plus */
	private ImagePlus imRef = null;

	/** working image width */
	private int width = 0;
	/** working image height */
	private int height = 0;
	/** working image depth */
	private int depth = 0;
	/** working image stack */
	private ImageStack inputImage = null;

	/** visit flags */
	private boolean [][][] visited = null;

	// Measures
	/** total number of end points voxels */
	private int totalNumberOfEndPoints = 0;
	/** total number of junctions voxels */
	private int totalNumberOfJunctionVoxels = 0;
	/** total number of slab voxels */
	private int totalNumberOfSlabs = 0;
	/** auxiliary variable to store the length of longest shortest path in a graph */
	private double shortestPath = 0;

	// Shortest path variables
	/** list of longest shortest paths from the skeletons in the image */
	private ArrayList< Double > shortestPathList;
	/** list containing longest shortest path points (one per tree) */
	private ArrayList< Point >[] shortestPathPoints = null; 
	/** shortest path x start position */
	private int spx = 0;
	/** shortest path y start position */
	private int spy = 0;
	/** shortest path z start position */
	private int spz = 0;
	/** shortest path start position array */
	private double[][] spStartPosition;
	/** shortest path output stack */
	private ImageStack shortPathImage = null;

	// Tree fields
	/** image stack containing all skeletons marked with their corresponding id */
	ImageStack labeledSkeletons = null;

	/** number of branches for every specific tree */
	private int[] numberOfBranches = null;
	/** number of end points voxels of every tree */
	private int[] numberOfEndPoints = null;
	/** number of junctions voxels of every tree*/
	private int[] numberOfJunctionVoxels = null;
	/** number of slab voxels of every specific tree */
	private int[] numberOfSlabs = null;
	/** number of junctions of every specific tree*/
	private int[] numberOfJunctions = null;
	/** number of triple points in every tree */
	private int[] numberOfTriplePoints = null;
	/** number of quadruple points in every tree */
	private int[] numberOfQuadruplePoints = null;
	/** list of end points in every tree */
	private ArrayList <Point> endPointsTree [] = null;
	/** list of junction voxels in every tree */
	private ArrayList <Point> junctionVoxelTree [] = null;
	/** list of special slab coordinates where circular tree starts */
	private ArrayList <Point> startingSlabTree [] = null;

	/** average branch length */
	private double[] averageBranchLength = null;

	/** maximum branch length */
	private double[] maximumBranchLength = null;

	/** list of end point coordinates in the entire image */
	private ArrayList <Point> listOfEndPoints = null;
	/** list of junction coordinates in the entire image */
	private ArrayList <Point> listOfJunctionVoxels = null;
	/** list of slab coordinates in the entire image */
	private ArrayList <Point> listOfSlabVoxels = null;
	/** list of slab coordinates in the entire image */
	private ArrayList <Point> listOfStartingSlabVoxels = null;

	/** list of groups of junction voxels that belong to the same tree junction (in every tree) */
	private ArrayList < ArrayList <Point> > listOfSingleJunctions[] = null;
	/** array of junction vertex per tree */
	private Vertex[][] junctionVertex = null;

	/** stack image containing the corresponding skeleton tags (end point, junction or slab) */
	private ImageStack taggedImage = null;

	/** auxiliary temporary point */
	private Point auxPoint = null;

	/** number of trees (skeletons) in the image */
	private int numOfTrees = 0;

	/** pruning option */
	private boolean bPruneCycles = true;

	 /** dead-end pruning option */
	public static boolean pruneEnds = DEFAULT_PRUNE_ENDS;

	/** protective-ROI option (branches inside ROI are spared from pruning) */
	public static boolean protectRoi = DEFAULT_PROTECT_ROI;

	/** calculate largest shortest path option */
	public static boolean calculateShortestPath = DEFAULT_CALCULATE_SHORTEST_PATH;

	/** array of graphs (one per tree) */
	private Graph[] graph = null;

	/** auxiliary list of slabs */
	private ArrayList<Point> slabList = null;
	/** auxiliary final vertex */
	private Vertex auxFinalVertex = null;

	/** prune cycle options */
	public static final String[] pruneCyclesModes = {"none", 
													 "shortest branch", 
													 "lowest intensity voxel", 
													 "lowest intensity branch"};
	/** no pruning mode index */
	public static final int NONE = 0;
	/** shortest branch pruning mode index */
	public static final int SHORTEST_BRANCH = 1;
	/** lowest pixel intensity pruning mode index */
	public static final int LOWEST_INTENSITY_VOXEL = 2;
	/** lowest intensity branch pruning mode index */
	public static final int LOWEST_INTENSITY_BRANCH = 3;

	/** original grayscale image (for lowest pixel intensity pruning mode) */
	private ImageStack originalImage = null;

	/** prune cycle options index */
	public static int pruneIndex = DEFAULT_PRUNE_MODE_INDEX;

	/** x- neighborhood offset */
	private int x_offset = 1;
	/** y- neighborhood offset */
	private int y_offset = 1;
	/** z- neighborhood offset */
	private int z_offset = 1;

	/** boolean flag to display extra information in result tables */
	public static boolean verbose = DEFAULT_VERBOSE;

	/** silent run flag, to distinguish between GUI and plugin calls */
	protected boolean silent = false;

	/** debugging flag */
	private static final boolean debug = false;

	/** flag to output the labeled skeletons in a new image */
	public static boolean displaySkeletons = DEFAULT_DISPLAY_SKELETONS;

	/* -----------------------------------------------------------------------*/
	/**
	 * This method is called once when the filter is loaded.
	 * 
	 * @param arg argument specified for this plugin
	 * @param imp currently active image
	 * @return flag word that specifies the filters capabilities
	 */
	public int setup(String arg, ImagePlus imp) 
	{
		this.imRef = imp;

		if (arg.equals("about")) 
		{
			showAbout();
			return DONE;
		}

		return DOES_8G; 
	} // end method setup 

	/* -----------------------------------------------------------------------*/
	/**
	 * Process the image: tag skeleton and show results.
	 * 
	 * @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
	 */
	public void run(ImageProcessor ip) 
	{
		loadDialogSettings();
		createSettingsDialog();
		settingsDialog.showDialog();
		if (settingsDialog.wasCanceled()) {
			return;
		}
		setSettingsFromDialog();
		saveDialogSettings();

		// pre-checking if another image is needed and also setting bPruneCycles
		ImagePlus origIP = null;
		switch(pruneIndex)
		{
			// No pruning
			case AnalyzeSkeleton_.NONE:
				this.bPruneCycles = false;
				break;
			// Pruning cycles by shortest branch
			case AnalyzeSkeleton_.SHORTEST_BRANCH:
				this.bPruneCycles = true;
				break;
			// Pruning cycles by lowest pixel intensity
			case AnalyzeSkeleton_.LOWEST_INTENSITY_VOXEL:
			case AnalyzeSkeleton_.LOWEST_INTENSITY_BRANCH:
				// Select original image between the open images
				int[] ids = WindowManager.getIDList();
				if ( ids == null || ids.length < 1 )
				{
					IJ.showMessage( "You should have at least one image open." );
					return;
				}

				String[] titles = new String[ ids.length ];
				for ( int i = 0; i < ids.length; ++i )
				{
					titles[ i ] = ( WindowManager.getImage( ids[ i ] ) ).getTitle();
				}

				final GenericDialog gd2 = new GenericDialog( "Image selection" );

				gd2.addMessage( "Select original grayscale image:" );
				final String current = WindowManager.getCurrentImage().getTitle();
				gd2.addChoice( "original_image", titles, current );

				gd2.showDialog();

				if (gd2.wasCanceled()) 
					return;

				// Get original stack
				origIP = WindowManager.getImage( ids[ gd2.getNextChoiceIndex() ] );

				this.bPruneCycles = true;
				break;
			default:
		}

		// now we have all the information that's needed for running the plugin
		// as if it was called from somewhere else
		run(pruneIndex, pruneEnds, calculateShortestPath, origIP, false,
				verbose, (protectRoi) ? this.imRef.getRoi() : null);

		if(debug)
			IJ.log("num of skeletons = " + this.numOfTrees);

		// Show labeled skeletons
		if( AnalyzeSkeleton_.displaySkeletons )
		{
			ImagePlus labeledSkeletons = 
					new ImagePlus( this.imRef.getShortTitle() 
							+ "-labeled-skeletons", this.labeledSkeletons.duplicate() );
			IJ.run( labeledSkeletons, "Fire", null );
			labeledSkeletons.show();
		}

		// Show results table
		showResults();

	} // end run method

	/** Disables dialog components that are irrelevant to GUI-based analysis. */
	public boolean dialogItemChanged(GenericDialog gd, AWTEvent e)
	{
		if ( this.imRef.getRoi() == null && null != gd
				&& null != gd.getCheckboxes() )
		{
			Checkbox roiOption = (Checkbox)gd.getCheckboxes().elementAt(1);
			roiOption.setEnabled(false);
			if (Recorder.record) roiOption.setState(false);
		}
		return true;
	}

	/**
	 * This method is intended for non-interactively using this plugin.
	 * <p>
	 * @param pruneIndex The pruneIndex, as asked by the initial gui dialog.
	 * @param pruneEnds flag to prune end-point-ending branches
	 * @param shortPath flag to calculate the longest shortest path
	 * @param origIP original grayscale input image (for lowest pixel intensity pruning mode)
	 * @param silent 
	 * @param verbose flag to display running information 
	 */
	public SkeletonResult run(
			int pruneIndex,
			boolean pruneEnds,
			boolean shortPath,
			ImagePlus origIP,
			boolean silent,
			boolean verbose)
	{
		return run(pruneIndex, pruneEnds, shortPath, origIP, silent, verbose, null);
	}

	/**
	 * This method is intended for non-interactively using this plugin.
	 * <p>
	 * @param pruneIndex The pruneIndex, as asked by the initial gui dialog.
	 * @param pruneEnds flag to prune end-point-ending branches
	 * @param shortPath flag to calculate the longest shortest path
	 * @param origIP original grayscale input image (for lowest pixel intensity pruning mode)
	 * @param silent
	 * @param verbose flag to display running information
	 * @param roi points inside this region are spared from elimination when
	 *            pruning end branches. ROI can be associated to a single image
	 *            in the stack or all images as per {@link ij.gui.Roi#getPosition
	 *            ij.gui.Roi.getPosition()}
	 */
	public SkeletonResult run(
			int pruneIndex,
			boolean pruneEnds,
			boolean shortPath,
			ImagePlus origIP,
			boolean silent,
			boolean verbose,
			Roi roi)
	{
		AnalyzeSkeleton_.pruneIndex = pruneIndex;
		this.silent = silent;
		AnalyzeSkeleton_.pruneEnds = pruneEnds;
		AnalyzeSkeleton_.calculateShortestPath = shortPath;
		AnalyzeSkeleton_.verbose = verbose;

		switch(pruneIndex)
		{
			// No pruning
			case AnalyzeSkeleton_.NONE:
				this.bPruneCycles = false;
				break;
			// Pruning cycles by shortest branch
			case AnalyzeSkeleton_.SHORTEST_BRANCH:
				this.bPruneCycles = true;
				break;
			// Pruning cycles by lowest pixel intensity
			case AnalyzeSkeleton_.LOWEST_INTENSITY_VOXEL:
			case AnalyzeSkeleton_.LOWEST_INTENSITY_BRANCH:
				// calculate neighborhood size given the calibration
				calculateNeighborhoodOffsets(origIP.getCalibration());
				this.originalImage = origIP.getStack();
				this.bPruneCycles = true;
				break;
			default:
		}

		this.width = this.imRef.getWidth();
		this.height = this.imRef.getHeight();
		this.depth = this.imRef.getStackSize();
		this.inputImage = this.imRef.getStack();

		// initialize visit flags
		resetVisited();

		// Tag skeleton, differentiate trees and visit them
		processSkeleton(this.inputImage);

		 // prune ends
		if (pruneEnds) 
		{
			pruneEndBranches(this.inputImage, this.taggedImage, roi);
		}

		// Prune cycles if necessary
		if(bPruneCycles)
		{
			if(pruneCycles(this.inputImage, this.originalImage, AnalyzeSkeleton_.pruneIndex))
			{
				// initialize visit flags
				resetVisited();
				// Recalculate analysis over the new image
				bPruneCycles = false;
				processSkeleton(this.inputImage);
			}
		}

		// Calculate triple points (junctions with exactly 3 branches)
		calculateTripleAndQuadruplePoints();

		if(shortPath)
		{
			if(debug)
				IJ.log("Calculating longest shortest paths...");

			// Copy input image
			this.shortPathImage = new ImageStack(this.width, this.height, this.inputImage.getColorModel());
			for(int i=1; i<=this.inputImage.getSize(); i++)
				shortPathImage.addSlice(this.inputImage.getSliceLabel(i), this.inputImage.getProcessor(i).duplicate());

			shortestPathList = new ArrayList < Double >();

			this.shortestPathPoints = new ArrayList [ this.numOfTrees ];

			// Visit skeleton and measure distances.
			// and apply Warshall algorithm
			spStartPosition = new double[this.numOfTrees][3];
			for(int i = 0; i < this.numOfTrees; i++)
			{
				shortestPathPoints[ i ] = new ArrayList<Point>();
				// Warshall algorithm including tag positions
				this.shortestPath = warshallAlgorithm(this.graph[i], shortestPathPoints[ i ]);
				shortestPathList.add(this.shortestPath);
				spStartPosition[i][0] = spx * this.imRef.getCalibration().pixelWidth;
				spStartPosition[i][1] = spy * this.imRef.getCalibration().pixelHeight;
				spStartPosition[i][2] = spz * this.imRef.getCalibration().pixelDepth;
			}

			if (!silent) {
				// Display short paths in a new stack
				ImagePlus shortIP = new ImagePlus("Longest shortest paths", shortPathImage);
				shortIP.show();

				// Set same calibration as the input image
				shortIP.setCalibration(this.imRef.getCalibration());

				// We apply the Fire LUT and reset the min and max to be between 0-255.
				IJ.run(shortIP, "Fire", null);

				//IJ.resetMinAndMax();
				shortIP.resetDisplayRange();
				shortIP.updateAndDraw();
			}

		}


		// Return the analysis results
		return assembleResults();
	}

	/**
	 * This method is intended for non-interactively using this plugin.
	 * <p>
	 * @param pruneIndex The pruneIndex, as asked by the initial gui dialog.
	 * @param thresholdLength maximum length of the branches to prune (all branches below that value are removed)
	 * @param shortPath flag to calculate the longest shortest path
	 * @param origIP original input image
	 * @param silent 
	 * @param verbose flag to display running information 
	 */
	public SkeletonResult run(
			int pruneIndex,
			double thresholdLength,
			boolean shortPath,
			ImagePlus origIP,
			boolean silent,
			boolean verbose)
	{
		AnalyzeSkeleton_.pruneIndex = pruneIndex;
		this.silent = silent;
		AnalyzeSkeleton_.pruneEnds = true;
		AnalyzeSkeleton_.calculateShortestPath = shortPath;
		AnalyzeSkeleton_.verbose = verbose;

		switch(pruneIndex)
		{
			// No pruning
			case AnalyzeSkeleton_.NONE:
				this.bPruneCycles = false;
				break;
			// Pruning cycles by shortest branch
			case AnalyzeSkeleton_.SHORTEST_BRANCH:
				this.bPruneCycles = true;
				break;
			// Pruning cycles by lowest pixel intensity
			case AnalyzeSkeleton_.LOWEST_INTENSITY_VOXEL:
			case AnalyzeSkeleton_.LOWEST_INTENSITY_BRANCH:
				// calculate neighborhood size given the calibration
				calculateNeighborhoodOffsets(origIP.getCalibration());
				this.originalImage = origIP.getStack();
				this.bPruneCycles = true;
				break;
			default:
		}

		this.width = this.imRef.getWidth();
		this.height = this.imRef.getHeight();
		this.depth = this.imRef.getStackSize();
		this.inputImage = this.imRef.getStack();

		// initialize visit flags
		resetVisited();

		// Tag skeleton, differentiate trees and visit them
		processSkeleton(this.inputImage);

		 // prune ends
		pruneEndBranches(this.inputImage, this.taggedImage, thresholdLength);


		// Prune cycles if necessary
		if(bPruneCycles)
		{
			if(pruneCycles(this.inputImage, this.originalImage, AnalyzeSkeleton_.pruneIndex))
			{
				// initialize visit flags
				resetVisited();
				// Recalculate analysis over the new image
				bPruneCycles = false;
				processSkeleton(this.inputImage);
			}
		}

		// Calculate triple points (junctions with exactly 3 branches)
		calculateTripleAndQuadruplePoints();

		if(shortPath)
		{
			if(debug)
				IJ.log("Calculating longest shortest paths...");

			// Copy input image
			this.shortPathImage = new ImageStack(this.width, this.height, this.inputImage.getColorModel());
			for(int i=1; i<=this.inputImage.getSize(); i++)
				shortPathImage.addSlice(this.inputImage.getSliceLabel(i), this.inputImage.getProcessor(i).duplicate());

			shortestPathList = new ArrayList < Double >();

			this.shortestPathPoints = new ArrayList [ this.numOfTrees ];

			// Visit skeleton and measure distances.
			// and apply Warshall algorithm
			spStartPosition = new double[this.numOfTrees][3];
			for(int i = 0; i < this.numOfTrees; i++)
			{
				shortestPathPoints[ i ] = new ArrayList<Point>();
				// Warshall algorithm including tag positions
				this.shortestPath = warshallAlgorithm(this.graph[i], shortestPathPoints[ i ]);
				shortestPathList.add(this.shortestPath);
				spStartPosition[i][0] = spx * this.imRef.getCalibration().pixelWidth;
				spStartPosition[i][1] = spy * this.imRef.getCalibration().pixelHeight;
				spStartPosition[i][2] = spz * this.imRef.getCalibration().pixelDepth;
			}

			if (!silent) {
				// Display short paths in a new stack
				ImagePlus shortIP = new ImagePlus("Longest shortest paths", shortPathImage.duplicate());
				shortIP.show();

				// Set same calibration as the input image
				shortIP.setCalibration(this.imRef.getCalibration());

				// We apply the Fire LUT and reset the min and max to be between 0-255.
				IJ.run(shortIP, "Fire", null);

				//IJ.resetMinAndMax();
				shortIP.resetDisplayRange();
				shortIP.updateAndDraw();
			}

		}


		// Return the analysis results
		return assembleResults();
	}

	/**
	 * Get the graphs of the current skeletons
	 * @return array of graphs (one per tree/skeleton)
	 */
	public Graph[] getGraphs()
	{
		return graph;
	}

	/**
	 * Get the list of points (including junctions and end points) of the 
	 * largest shortest paths in the skeleton image (one per tree).
	 * @return array with the lists of points of the shortest paths
	 */
	public ArrayList<Point>[] getShortestPathPoints()
	{
		return this.shortestPathPoints;
	}
		/**
	 * Get the list of endpoints in the skeleton image per tree.
	 * @return array with the lists of endpoints
	 */
	public ArrayList<Point>[] getEndPointsTree()
	{
		return this.endPointsTree;
	}
	/**
	 * Get the list of endpoints in the skeleton image.
	 * @return array with the endpoints
	 */
	public ArrayList<Point> getEndPoints()
	{
		return this.listOfEndPoints;
	}


	/**
	 * A simpler standalone running method, for analysis without pruning
	 * or showing images.
	 * <p>
	 * This one just calls run(AnalyzeSkeleton_.NONE, false, null, true, false)
	 */
	public SkeletonResult run()
	{
		return run(NONE, false, false, null, true, false);
	}

	/**
	 * Prune end branches outside the specified ROI
	 *
	 * @param stack input skeleton image
	 * @param taggedImage tagged skeleton image
	 * @param roi 'protective' ROI: points inside this region are spared from pruning
	 *
	 */
	private void pruneEndBranches(ImageStack stack, ImageStack taggedImage, Roi roi)
	{
		if(debug)
			IJ.log("Pruning end-point branches...");
		for (int t = 0; t < this.numOfTrees; t++)
		{
			if(debug)
				IJ.log("Pruning tree #" + t);

			Graph g = graph[t];
			ArrayList<Vertex> vertices = g.getVertices();
			ListIterator<Vertex> vit = vertices.listIterator();

			if(debug)
				IJ.log("Initial number of vertices: " + graph[t].getVertices().size());

			while (vit.hasNext())
			{
				Vertex v = vit.next();
				// Check if the vertex is an end point
				if (v.getBranches().size() == 1 && isEndPoint( v.getPoints().get( 0 ), roi) )
				{
					if(debug)
						IJ.log("Pruning branch starting at " + v.getPoints().get(0));

					// Remove end point voxels
					ArrayList<Point> points = v.getPoints();
					final int nPoints = points.size();

					for (int i = 0; i < nPoints; i++)
					{
						Point p = points.get(i);
						setPixel(stack, p.x, p.y, p.z, (byte) 0);
						setPixel(taggedImage, p.x, p.y, p.z, (byte) 0);
						this.numberOfEndPoints[t]--;
						this.totalNumberOfEndPoints--;
						Iterator<Point> pit = this.listOfEndPoints.listIterator();
						while (pit.hasNext()){
							Point ep = pit.next();
							if (ep.equals(p)){
								pit.remove();
								break;
							}
						}
					}

					// Remove branch voxels
					Edge branch = v.getBranches().get(0);
					points = branch.getSlabs();
					final int nSlabs = points.size();
					for (int i = 0; i < nSlabs; i++)
					{
						Point p = points.get(i);
						setPixel(stack, p.x, p.y, p.z, (byte) 0);
						setPixel(taggedImage, p.x, p.y, p.z, (byte) 0);
						this.numberOfSlabs[t]--;
						this.totalNumberOfSlabs--;
						Iterator<Point> pit = this.listOfSlabVoxels.listIterator();
						while (pit.hasNext())
						{
							Point ep = pit.next();
							if (ep.equals(p)){
								pit.remove();
								break;
							}
						}
					}

					// remove the Edge from the Graph
					ArrayList<Edge> gEdges = graph[t].getEdges();
					Iterator<Edge> git = gEdges.listIterator();
					while (git.hasNext())
					{
						Edge e = git.next();
						if (e.equals(branch))
						{
							git.remove();
							break;
						}
					}

					// remove the Edge from the opposite Vertex
					Vertex opp = branch.getOppositeVertex(v);
					ArrayList<Edge> oppBranches = opp.getBranches();
					Iterator<Edge> oppIt = oppBranches.listIterator();
					while (oppIt.hasNext())
					{
						Edge oppBranch = oppIt.next();
						if (oppBranch.equals(branch))
						{
							oppIt.remove();
							break;
						}
					}

					// remove the Edge from the Vertex
					v.getBranches().remove(0);

					// remove the Vertex from the Graph
					vit.remove();
				}
			}

			if(debug)
				IJ.log("Final number of vertices: " + graph[t].getVertices().size());
		}

		return;
	}

	/**
	 * Prune end branches of a specific length
	 *
	 * @param stack input skeleton image
	 * @param taggedImage tagged skeleton image
	 * @param length limit length to prune the branches (in calibrated units)
	 *
	 */
	private void pruneEndBranches(
			ImageStack stack, 
			ImageStack taggedImage, 
			double length) 
	{
		if(debug)
			IJ.log("Pruning end-point branches...");
		for (int t = 0; t < this.numOfTrees; t++)
		{
			if(debug)
				IJ.log("Pruning tree #" + t);

			Graph g = graph[t];
			ArrayList<Vertex> vertices = g.getVertices();
			ListIterator<Vertex> vit = vertices.listIterator();

			if(debug)
				IJ.log("Initial number of vertices: " + graph[t].getVertices().size());

			while (vit.hasNext())
			{
				Vertex v = vit.next();
				if (v.getBranches().size() == 1 && v.getBranches().get(0).getLength() <= length)
				{
					if(debug)
						IJ.log("Pruning branch starting at " + v.getPoints().get(0));
					// Remove end point voxels
					ArrayList<Point> points = v.getPoints();
					final int nPoints = points.size();

					for (int i = 0; i < nPoints; i++)
					{
						Point p = points.get(i);
						setPixel(stack, p.x, p.y, p.z, (byte) 0);
						setPixel(taggedImage, p.x, p.y, p.z, (byte) 0);
						this.numberOfEndPoints[t]--;
						this.totalNumberOfEndPoints--;
						Iterator<Point> pit = this.listOfEndPoints.listIterator();
						while (pit.hasNext()){
							Point ep = pit.next();
							if (ep.equals(p)){
								pit.remove();
								break;
							}
						}
					}

					// Remove branch voxels
					Edge branch = v.getBranches().get(0);
					points = branch.getSlabs();
					final int nSlabs = points.size();
					for (int i = 0; i < nSlabs; i++)
					{
						Point p = points.get(i);
						setPixel(stack, p.x, p.y, p.z, (byte) 0);
						setPixel(taggedImage, p.x, p.y, p.z, (byte) 0);
						this.numberOfSlabs[t]--;
						this.totalNumberOfSlabs--;
						Iterator<Point> pit = this.listOfSlabVoxels.listIterator();
						while (pit.hasNext())
						{
							Point ep = pit.next();
							if (ep.equals(p)){
								pit.remove();
								break;
							}
						}
					}

					// remove the Edge from the Graph
					ArrayList<Edge> gEdges = graph[t].getEdges();
					Iterator<Edge> git = gEdges.listIterator();
					while (git.hasNext())
					{
						Edge e = git.next();
						if (e.equals(branch))
						{
							git.remove();
							break;
						}
					}

					// remove the Edge from the opposite Vertex
					Vertex opp = branch.getOppositeVertex(v);
					ArrayList<Edge> oppBranches = opp.getBranches();
					Iterator<Edge> oppIt = oppBranches.listIterator();
					while (oppIt.hasNext())
					{
						Edge oppBranch = oppIt.next();
						if (oppBranch.equals(branch))
						{
							oppIt.remove();
							break;
						}
					}

					// remove the Edge from the Vertex
					v.getBranches().remove(0);

					// remove the Vertex from the Graph
					vit.remove();
				}
			}

			if(debug)
				IJ.log("Final number of vertices: " + graph[t].getVertices().size());
		}

		return;
	}


	// ---------------------------------------------------------------------------
	/**
	 * Calculate the neighborhood size based on the calibration of the image.
	 * @param calibration image calibration
	 */
	private void calculateNeighborhoodOffsets(Calibration calibration) 
	{
		double max = calibration.pixelDepth;
		if(calibration.pixelHeight > max)
			max = calibration.pixelHeight;
		if(calibration.pixelWidth > max)
			max = calibration.pixelWidth;

		this.x_offset = ((int) Math.round(max / calibration.pixelWidth) > 1) ? 
							(int) Math.round(max / calibration.pixelWidth) : 1;
		this.y_offset = ((int) Math.round(max / calibration.pixelHeight) > 1) ? 
							(int) Math.round(max / calibration.pixelHeight) : 1;
		this.z_offset = ((int) Math.round(max / calibration.pixelDepth) > 1) ? 
							(int) Math.round(max / calibration.pixelDepth) : 1;

		if(debug)
		{
			IJ.log("x_offset = " + this.x_offset);
			IJ.log("y_offset = " + this.y_offset);
			IJ.log("z_offset = " + this.z_offset);
		}

	}// end method calculateNeighborhoodOffsets

	// ---------------------------------------------------------------------------
	/**
	 * Process skeleton: tag image, mark trees and visit.
	 * 
	 * @param inputImage2 input skeleton image to process
	 */
	public void processSkeleton(ImageStack inputImage2) 
	{
		// Initialize  global lists of points
		this.listOfEndPoints = new ArrayList<Point>();
		this.listOfJunctionVoxels = new ArrayList<Point>();
		this.listOfSlabVoxels = new ArrayList<Point>();
		this.listOfStartingSlabVoxels = new ArrayList<Point>();
		this.totalNumberOfEndPoints = 0;
		this.totalNumberOfJunctionVoxels = 0;
		this.totalNumberOfSlabs = 0;


		// Prepare data: classify voxels and tag them.
		this.taggedImage = tagImage(inputImage2);

		// Show tags image.
		if(!bPruneCycles && !silent)
		{
			displayTagImage(taggedImage);
		}

		// Mark trees
		labeledSkeletons = markTrees(taggedImage);

		if(this.numOfTrees == 0)
			return;

		// Ask memory for every tree
		initializeTrees();

		// Divide groups of end-points and junction voxels
		if(this.numOfTrees > 1)
			divideVoxelsByTrees( labeledSkeletons );
		if(this.numOfTrees == 1)
		{
			if(debug)
				IJ.log("list of end points size = " + this.listOfEndPoints.size());
			this.endPointsTree[0] = this.listOfEndPoints;
			this.numberOfEndPoints[0] = this.listOfEndPoints.size();
			this.junctionVoxelTree[0] = this.listOfJunctionVoxels;
			this.numberOfJunctionVoxels[0] = this.listOfJunctionVoxels.size();
			this.startingSlabTree[0] = this.listOfStartingSlabVoxels;
		}

		// Calculate number of junctions (skipping neighbor junction voxels)
		groupJunctions( labeledSkeletons );

		// Mark all unvisited
		resetVisited();

		// Visit skeleton and measure distances.
		for(int i = 0; i < this.numOfTrees; i++)
			visitSkeleton(taggedImage, labeledSkeletons, i+1);

	} // end method processSkeleton

	// -----------------------------------------------------------------------
	/**
	 * Prune cycles from tagged image and update it.
	 * 
	 * @param inputImage input skeleton image
	 * @param originalImage original gray-scale image
	 * @param pruningMode (SHORTEST_BRANCH, LOWEST_INTENSITY_VOXEL, LOWEST_INTENSITY_BRANCH)
	 * @return true if the input image was pruned or false if there were no cycles
	 */
	private boolean pruneCycles(
			ImageStack inputImage, 
			final ImageStack originalImage, 
			final int pruningMode) 
	{
		boolean pruned = false;

		for(int iTree = 0 ; iTree < this.numOfTrees; iTree ++)
		{
			// For circular trees we just remove one slab
			if(this.startingSlabTree[iTree].size() == 1)
			{
				setPixel(inputImage, this.startingSlabTree[iTree].get(0),(byte) 0);
				pruned = true;
			}
			else // For the rest, we do depth-first search to detect the cycles
			{
				// DFS
				ArrayList <Edge> backEdges = this.graph[iTree].depthFirstSearch();

				if(debug)
				{
					IJ.log( " --------------------------- ");
					final String[] s = new String[]{"UNDEFINED", "TREE" , "BACK"};
					for(final Edge e : this.graph[iTree].getEdges())
					{
						IJ.log(" edge " + e.getV1().getPoints().get(0) + " - " + e.getV2().getPoints().get(0) + " : " + s[e.getType()+1]);
					}
				}

				// If DFS returned backEdges, we need to delete the loops
				if(backEdges.size() > 0)
				{
					// Find all edges of each loop (backtracking the predecessors)
					for(final Edge e : backEdges)
					{
						ArrayList<Edge> loopEdges = new ArrayList<Edge>();
						loopEdges.add(e);

						Edge minEdge = e;

						// backtracking (starting at the vertex with higher order index
						final Vertex finalLoopVertex = e.getV1().getVisitOrder() < e.getV2().getVisitOrder() ? e.getV1() : e.getV2();

						Vertex backtrackVertex = e.getV1().getVisitOrder() < e.getV2().getVisitOrder() ? e.getV2() : e.getV1();

						// backtrack until reaching final loop vertex
						while(!finalLoopVertex.equals(backtrackVertex))
						{
							// Extract predecessor
							final Edge pre = backtrackVertex.getPredecessor();
							// Update shortest loop edge if necessary
							if(pruningMode == AnalyzeSkeleton_.SHORTEST_BRANCH && 
									pre.getSlabs().size() < minEdge.getSlabs().size())
								minEdge = pre;
							// Add to loop edge list
							loopEdges.add(pre);
							// Extract predecessor
							backtrackVertex = pre.getV1().equals(backtrackVertex) ? pre.getV2() : pre.getV1(); 
						}

						// Prune cycle
						if(pruningMode == AnalyzeSkeleton_.SHORTEST_BRANCH)
						{
							// Remove middle slab from the shortest loop edge
							Point removeCoords = null;
							if(minEdge.getSlabs().size() > 0)
								removeCoords = minEdge.getSlabs().get(minEdge.getSlabs().size()/2);
							else 
								removeCoords = minEdge.getV1().getPoints().get(0);
							setPixel(inputImage, removeCoords,(byte) 0);
						}
						else if (pruningMode == AnalyzeSkeleton_.LOWEST_INTENSITY_VOXEL)
						{
							removeLowestIntensityVoxel(loopEdges, inputImage, originalImage);
						}
						else if(pruningMode == AnalyzeSkeleton_.LOWEST_INTENSITY_BRANCH)
						{
							cutLowestIntensityBranch(loopEdges, inputImage, originalImage);
						}
					}// endfor backEdges

					pruned = true;
				}
			}
		}

		return pruned;
	}// end method pruneCycles

	// -----------------------------------------------------------------------
	/**
	 * Cut the a list of edges in the lowest pixel intensity voxel (calculated
	 * from the original -grayscale- image).
	 * 
	 * @param loopEdges list of edges to be analyzed
	 * @param inputImage2 input skeleton image
	 * @param originalGrayImage original gray image
	 */
	private void removeLowestIntensityVoxel(
			final ArrayList<Edge> loopEdges,
			ImageStack inputImage2,
			ImageStack originalGrayImage) 
	{
		Point lowestIntensityVoxel = null;

		double lowestIntensityValue = Double.MAX_VALUE;

		for(final Edge e : loopEdges)
		{
			for(final Point p : e.getSlabs())
			{
				final double avg = getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				if(avg < lowestIntensityValue)
				{
					lowestIntensityValue = avg;
					lowestIntensityVoxel = p;
				}
			}
			// Check vertices
			/*
			for(final Point p : e.getV1().getPoints())
			{
				final double avg = getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				if(avg < lowestIntensityValue)
				{
					lowestIntensityValue = avg;
					lowestIntensityVoxel = p;
				}
			}
			for(final Point p : e.getV2().getPoints())
			{
				final double avg = getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				if(avg < lowestIntensityValue)
				{
					lowestIntensityValue = avg;
					lowestIntensityVoxel = p;
				}
			}
			*/
		}

		// Cut loop in the lowest intensity pixel value position
		if(debug)
			IJ.log("Cut loop at coordinates: " + lowestIntensityVoxel);
		setPixel(inputImage2, lowestIntensityVoxel,(byte) 0);
	}//end method removeLowestIntensityVoxel

	// -----------------------------------------------------------------------
	/**
	 * Cut the a list of edges in the lowest pixel intensity branch.
	 * 
	 * @param loopEdges list of edges to be analyzed
	 * @param inputImage2 input skeleton image
	 */
	private void cutLowestIntensityBranch(
			final ArrayList<Edge> loopEdges,
			ImageStack inputImage2,
			ImageStack originalGrayImage) 
	{
		Edge lowestIntensityEdge = null;

		double lowestIntensityValue = Double.MAX_VALUE;

		Point cutPoint = null;

		for(final Edge e : loopEdges)
		{
			// Calculate average intensity of the edge neighborhood
			double min_val = Double.MAX_VALUE;
			Point darkestPoint = null;

			double edgeIntensity = 0;
			double n_vox = 0;

			// Check slab points
			for(final Point p : e.getSlabs())
			{
				final double avg = getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				// Keep track of the darkest slab point of the edge
				if(avg < min_val)
				{
					min_val = avg;
					darkestPoint = p;
				}
				edgeIntensity += avg;
				n_vox++;
			}
			// Check vertices
			for(final Point p : e.getV1().getPoints())
			{
				edgeIntensity += getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				n_vox++;
			}
			for(final Point p : e.getV2().getPoints())
			{
				edgeIntensity += getAverageNeighborhoodValue(originalGrayImage, p,
						this.x_offset, this.y_offset, this.z_offset);
				n_vox++;
			}

			if(n_vox != 0)
				edgeIntensity /= n_vox;
			if(debug)
			{
				IJ.log("Loop edge between " + e.getV1().getPoints().get(0) + " and " + e.getV2().getPoints().get(0) + ":");
				IJ.log("avg edge intensity = " + edgeIntensity + " darkest slab point = " + darkestPoint.toString());
			}
			// Keep track of the lowest intensity edge
			if(edgeIntensity < lowestIntensityValue)
			{
				lowestIntensityEdge = e;
				lowestIntensityValue = edgeIntensity;
				cutPoint = darkestPoint;
			}
		}

		// Cut loop in the lowest intensity branch medium position
		Point removeCoords = null;
		if (lowestIntensityEdge.getSlabs().size() > 0)
			removeCoords = cutPoint;
		else 
		{
			IJ.error("Lowest intensity branch without slabs?!: vertex " 
					+ lowestIntensityEdge.getV1().getPoints().get(0));
			removeCoords = lowestIntensityEdge.getV1().getPoints().get(0);
		}

		if(debug)
			IJ.log("Cut loop at coordinates: " + removeCoords);
		setPixel(inputImage2, removeCoords,(byte) 0);

	}// end method cutLowestIntensityBranch

	// -----------------------------------------------------------------------
	/**
	 * Display tag image on a new window.
	 * 
	 * @param taggedImage tag image to be displayed
	 */
	void displayTagImage(ImageStack taggedImage)
	{
		final int slices = imRef.getNSlices();
		final int frames = imRef.getNFrames();
		final int channels = imRef.getNChannels();

		ImagePlus tagIP = IJ.createHyperStack("Tagged skeleton", width, height, channels, slices, frames,
											  imRef.getBitDepth());
		tagIP.setStack(taggedImage.duplicate(), channels, slices, frames);
		tagIP.show();

		// Set same calibration as the input image
		tagIP.setCalibration(this.imRef.getCalibration());

		// We apply the Fire LUT and reset the min and max to be between 0-255.
		IJ.run(tagIP, "Fire", null);

		//IJ.resetMinAndMax();
		tagIP.resetDisplayRange();
		tagIP.updateAndDraw();
	} // end method displayTagImage

	// -----------------------------------------------------------------------
	/**
	 * Divide the end point, junction and special (starting) slab voxels in the 
	 * corresponding tree lists.
	 * 
	 *  @param treeIS tree image
	 */
	private void divideVoxelsByTrees(ImageStack treeIS) 
	{
		// Add end points to the corresponding tree
		for(int i = 0; i < this.totalNumberOfEndPoints; i++)
		{
			final Point p = this.listOfEndPoints.get(i);
			this.endPointsTree[ (int)getFloatPixel( treeIS, p ) - 1 ].add(p);
		}

		// Add junction voxels to the corresponding tree
		for(int i = 0; i < this.totalNumberOfJunctionVoxels; i++)
		{
			final Point p = this.listOfJunctionVoxels.get(i);
			this.junctionVoxelTree[ (int)getFloatPixel( treeIS, p ) - 1 ].add(p);;
		}

		// Add special slab voxels to the corresponding tree
		for(int i = 0; i < this.listOfStartingSlabVoxels.size(); i++)
		{
			final Point p = this.listOfStartingSlabVoxels.get(i);
			this.startingSlabTree[ (int)getFloatPixel( treeIS, p ) - 1 ].add(p);;
		}

		// Assign number of end points and junction voxels per tree
		for(int iTree = 0; iTree < this.numOfTrees; iTree++)
		{
			this.numberOfEndPoints[iTree] = this.endPointsTree[iTree].size();
			this.numberOfJunctionVoxels[iTree] = this.junctionVoxelTree[iTree].size();
		}

	} // end divideVoxelsByTrees

	// -----------------------------------------------------------------------
	/**
	 * Ask memory for trees.
	 */
	private void initializeTrees()
	{
		this.numberOfBranches = new int[this.numOfTrees];
		this.numberOfEndPoints = new int[this.numOfTrees];
		this.numberOfJunctionVoxels = new int[this.numOfTrees];
		this.numberOfJunctions = new int[this.numOfTrees];
		this.numberOfSlabs = new int[this.numOfTrees];
		this.numberOfTriplePoints = new int[this.numOfTrees];
		this.numberOfQuadruplePoints = new int[this.numOfTrees];
		this.averageBranchLength = new double[this.numOfTrees];
		this.maximumBranchLength = new double[this.numOfTrees];
		this.endPointsTree = new ArrayList[this.numOfTrees];
		this.junctionVoxelTree = new ArrayList[this.numOfTrees];
		this.startingSlabTree = new ArrayList[this.numOfTrees];
		this.listOfSingleJunctions = new ArrayList[this.numOfTrees];

		this.graph = new Graph[this.numOfTrees];

		for(int i = 0; i < this.numOfTrees; i++)
		{
			this.endPointsTree[i] = new ArrayList <Point>();
			this.junctionVoxelTree[i] = new ArrayList <Point>();
			this.startingSlabTree[i] = new ArrayList <Point>();
			this.listOfSingleJunctions[i] = new ArrayList < ArrayList <Point> > ();
		}
		this.junctionVertex  = new Vertex[this.numOfTrees][];
	}// end method initializeTrees

	// -----------------------------------------------------------------------
	/**
	 * Show results table.
	 */
	private void showResults() 
	{
		final ResultsTable rt = new ResultsTable();
		rt.showRowNumbers( true );

		final String[] head = {"Skeleton", "# Branches","# Junctions", "# End-point voxels",
						 "# Junction voxels","# Slab voxels","Average Branch Length", 
						 "# Triple points", "# Quadruple points", "Maximum Branch Length",
						 "Longest Shortest Path", "spx", "spy", "spz"};

		for(int i = 0 ; i < this.numOfTrees; i++)
		{
			rt.incrementCounter();

			rt.addValue(head[ 1], this.numberOfBranches[i]);        
			rt.addValue(head[ 2], this.numberOfJunctions[i]);
			rt.addValue(head[ 3], this.numberOfEndPoints[i]);
			rt.addValue(head[ 4], this.numberOfJunctionVoxels[i]);
			rt.addValue(head[ 5], this.numberOfSlabs[i]);
			rt.addValue(head[ 6], this.averageBranchLength[i]);
			rt.addValue(head[ 7], this.numberOfTriplePoints[i]);
			rt.addValue(head[ 8], this.numberOfQuadruplePoints[i]);
			rt.addValue(head[ 9], this.maximumBranchLength[i]);
			if(null != this.shortestPathList)
			{
				rt.addValue(head[10],this.shortestPathList.get(i));
				rt.addValue(head[11],this.spStartPosition[i][0]);
				rt.addValue(head[12],this.spStartPosition[i][1]);
				rt.addValue(head[13],this.spStartPosition[i][2]);
			}

			if (0 == i % 100) 
				rt.show("Results");
		}
		rt.show("Results");

		// Extra information
		if(AnalyzeSkeleton_.verbose)
		{
			// New results table
			final ResultsTable extra_rt = new ResultsTable();
			rt.showRowNumbers( true );

			final String[] extra_head = {"Branch", "Skeleton ID", 
							"Branch length","V1 x", "V1 y",
							"V1 z","V2 x","V2 y", "V2 z", "Euclidean distance","running average length", "average intensity (inner 3rd)", "average intensity"};


			// Edge comparator (by branch length)
			Comparator<Edge> comp = new Comparator<Edge>(){
				public int compare(Edge o1, Edge o2)
				{
					final double diff = o1.getLength() - o2.getLength(); 
					if(diff < 0)
						return 1;
					else if(diff == 0)
						return 0;
					else
						return -1;
				}
				public boolean equals(Object o)
				{
					return false;
				}
			};
			// Display branch information for each tree
			for(int i = 0 ; i < this.numOfTrees; i++)
			{
				final ArrayList<Edge> listEdges = this.graph[i].getEdges();
				// Sort branches by length
				Collections.sort(listEdges, comp);
				for(final Edge e : listEdges)
				{
					extra_rt.incrementCounter();
					extra_rt.addValue(extra_head[1], i+1);
					extra_rt.addValue(extra_head[2], e.getLength());
					extra_rt.addValue(extra_head[3], e.getV1().getPoints().get(0).x * this.imRef.getCalibration().pixelWidth);
					extra_rt.addValue(extra_head[4], e.getV1().getPoints().get(0).y * this.imRef.getCalibration().pixelHeight);
					extra_rt.addValue(extra_head[5], e.getV1().getPoints().get(0).z * this.imRef.getCalibration().pixelDepth);
					extra_rt.addValue(extra_head[6], e.getV2().getPoints().get(0).x * this.imRef.getCalibration().pixelWidth);
					extra_rt.addValue(extra_head[7], e.getV2().getPoints().get(0).y * this.imRef.getCalibration().pixelHeight);
					extra_rt.addValue(extra_head[8], e.getV2().getPoints().get(0).z * this.imRef.getCalibration().pixelDepth);
					extra_rt.addValue(extra_head[9], this.calculateDistance(e.getV1().getPoints().get(0), e.getV2().getPoints().get(0)));
					extra_rt.addValue(extra_head[10], e.getLength_ra());
					extra_rt.addValue(extra_head[11], e.getColor3rd());
					extra_rt.addValue(extra_head[12], e.getColor());
				}

			}
			extra_rt.show("Branch information");
		}

	}// end method showResults

	/**
	 * Returns one of the two result images in an ImageStack object.
	 *
	 * @param longestShortestPath Get the tagged longest shortest paths instead of the standard tagged image
	 *
	 * @return The results image with a tagged skeleton 
	 */
	public ImageStack getResultImage(boolean longestShortestPath)
	{
		if (longestShortestPath) {
			return this.shortPathImage;
		}
		return this.taggedImage;
	}


	/**
	 * Returns the analysis results in a SkeletonResult object.
	 * <p>
	 *
	 * @return The results of the skeleton analysis.
	 */
	protected SkeletonResult assembleResults()
	{
		SkeletonResult result = new SkeletonResult(numOfTrees);
		result.setBranches(numberOfBranches);
		result.setJunctions(numberOfJunctions);
		result.setEndPoints(numberOfEndPoints);
		result.setJunctionVoxels(numberOfJunctionVoxels);
		result.setSlabs(numberOfSlabs);
		result.setAverageBranchLength(averageBranchLength);
		result.setTriples(numberOfTriplePoints);
		result.setQuadruples(numberOfQuadruplePoints);
		result.setMaximumBranchLength(maximumBranchLength);

		result.setListOfEndPoints(listOfEndPoints);
		result.setListOfJunctionVoxels(listOfJunctionVoxels);
		result.setListOfSlabVoxels(listOfSlabVoxels);
		result.setListOfStartingSlabVoxels(listOfStartingSlabVoxels);

		result.setShortestPathList(shortestPathList);
		result.setSpStartPosition(spStartPosition);

		result.setGraph(graph);

		result.calculateNumberOfVoxels();

		return result;
	}

	// -----------------------------------------------------------------------
	/**
	 * Visit skeleton starting at end-points, junctions and slab of circular 
	 * skeletons, and record measurements.
	 * 
	 * @param taggedImage tag skeleton image
	 * @param treeImage skeleton image with tree classification
	 * @param currentTree number of the tree to be visited
	 */
	private void visitSkeleton(ImageStack taggedImage, ImageStack treeImage, int currentTree) 
	{
		// tree index
		final int iTree = currentTree - 1;

		if(debug)
		{
			// Junction vertices in the tree
			IJ.log("this.junctionVertex["+(iTree)+"].length = " + this.junctionVertex[iTree].length);
			for(int i = 0; i < this.junctionVertex[iTree].length; i++)
			{
				IJ.log(" vertices points: " + this.junctionVertex[iTree][i]);
			}
		}

		// Create new graph
		this.graph[iTree] = new Graph();
		// Add all junction vertices
		for(int i = 0; i < this.junctionVertex[iTree].length; i++)
			this.graph[iTree].addVertex(this.junctionVertex[iTree][i]);

		if(debug)
			IJ.log(" Analyzing tree number " + currentTree);
		// length of branches
		double branchLength = 0;

		this.maximumBranchLength[iTree] = 0;
		this.numberOfSlabs[iTree] = 0;

		// Visit branches starting at end points
		for(int i = 0; i < this.numberOfEndPoints[iTree]; i++)
		{
			final Point endPointCoord = this.endPointsTree[iTree].get(i);

			if(debug)
				IJ.log("\n*** visit from end point: " + endPointCoord + " *** ");

			// Skip when visited
			if(isVisited(endPointCoord))
			{
				//if(this.initialPoint[iTree] == null)
				//	IJ.error("WEIRD:" + " (" + endPointCoord.x + ", " + endPointCoord.y + ", " + endPointCoord.z + ")");
				if(debug)
					IJ.log("visited = (" + endPointCoord.x + ", " + endPointCoord.y + ", " + endPointCoord.z + ")");
				continue;
			}

			// Initial vertex
			Vertex v1 = new Vertex();
			v1.addPoint(endPointCoord);
			this.graph[iTree].addVertex(v1);
			if(i == 0)
				this.graph[iTree].setRoot(v1);

			// slab list for the edge
			this.slabList = new ArrayList<Point>();

			// Otherwise, visit branch until next junction or end point.
			double[] properties = visitBranch(endPointCoord, iTree);
			double length = properties[0];
			double color3rd = properties[1];
			double color = properties[2];
			double length_ra = properties[3];


			// If length is 0, it means the tree is formed by only one voxel.
			if(length == 0)
			{
				// If there is an adjacent visited junction, count it
				// as a single voxel branch
				final Point aux = getVisitedJunctionNeighbor(endPointCoord, v1);
				if(null != aux)
				{
					this.auxFinalVertex = findPointVertex(this.junctionVertex[iTree], aux);
					length += calculateDistance(endPointCoord, aux);

					// Add the length to the first point of the vertex (to prevent later from having
					// euclidean distances larger than the actual distance)
					length += calculateDistance(auxFinalVertex.getPoints().get(0), endPointCoord);
					// Add branch to graph
					if(debug)
						IJ.log( "adding branch from " + v1.getPoints().get(0) + " to " + this.auxFinalVertex.getPoints().get(0) );
					this.graph[iTree].addVertex(this.auxFinalVertex);
					this.graph[iTree].addEdge(new Edge(v1, this.auxFinalVertex, this.slabList, length, color3rd, color, length_ra));
					// increase number of branches
					this.numberOfBranches[iTree]++;

					if(debug)
						IJ.log("increased number of branches, length = " + length);

					branchLength += length;
				}
				else
					if(debug)
						IJ.log("set initial point to final point");
				continue;
			}

			// If the final point is a slab, then we add the path to the
			// neighbor junction voxel not belonging to the initial vertex
			// (unless it is a self loop)
			if(isSlab(this.auxPoint))
			{
				final Point aux = this.auxPoint;
				//IJ.log("Looking for " + this.auxPoint + " in the list of vertices...");
				this.auxPoint = getVisitedJunctionNeighbor(this.auxPoint, v1);
				this.auxFinalVertex = findPointVertex(this.junctionVertex[iTree], this.auxPoint);
				if(this.auxPoint == null)
				{
					//IJ.log("Point "+ aux + " has not neighbor end junction! (inner loop)");
					// Inner loop
					this.auxFinalVertex = v1;
					this.auxPoint = aux;
				}
				length += calculateDistance(this.auxPoint, aux);

				// Add the length to the first point of the vertex (to prevent later from having
				// euclidean distances larger than the actual distance)
				length += calculateDistance(auxFinalVertex.getPoints().get(0), auxPoint);
			}

			// Add branch to graph
			if(debug)
				IJ.log("adding branch from " + v1.getPoints().get(0) + " to " + this.auxFinalVertex.getPoints().get(0) +  ", aux point = " + this.auxPoint);
			this.graph[iTree].addVertex(this.auxFinalVertex);
			this.graph[iTree].addEdge(new Edge(v1, this.auxFinalVertex, this.slabList, length, color3rd, color, length_ra));

			// increase number of branches
			this.numberOfBranches[iTree]++;

			if(debug)
				IJ.log("increased number of branches, length = " + length);

			branchLength += length;

			// update maximum branch length
			if(length > this.maximumBranchLength[iTree])
			{
				this.maximumBranchLength[iTree] = length;
			}
		}

		// If there is no end points, set the first junction as root.
		if(this.numberOfEndPoints[iTree] == 0 && this.junctionVoxelTree[iTree].size() > 0)
			this.graph[iTree].setRoot(this.junctionVertex[iTree][0]);

		if(debug)
			IJ.log( " --------------------------- ");

		// Now visit branches starting at junctions
		// 08/26/2009 Changed the loop to visit first the junction voxels that are
		//            forming a single junction.
		for(int i = 0; i < this.junctionVertex[iTree].length; i++)
		{
			for(int j = 0; j < this.junctionVertex[iTree][i].getPoints().size(); j++)
			{
				final Point junctionCoord = this.junctionVertex[iTree][i].getPoints().get(j);

				if(debug)
					IJ.log("\n*** visit from junction " + junctionCoord + " *** ");

				// Mark junction as visited
				setVisited(junctionCoord, true);

				Point nextPoint = getNextUnvisitedVoxel(junctionCoord);

				while(nextPoint != null)
				{
					// Do not count adjacent junctions
					if( !isJunction(nextPoint))
					{
						if (debug)
							IJ.log("visiting " + nextPoint);

						// Create graph edge
						this.slabList = new ArrayList<Point>();
						this.slabList.add(nextPoint);
						this.numberOfSlabs[iTree]++;

						// Calculate distance from junction to that point
						double length = calculateDistance(junctionCoord, nextPoint);

						// Visit branch
						this.auxPoint = null;

						double[] properties = visitBranch(nextPoint, iTree);
						length += properties[0];
						double color3rd = properties[1];
						double color = properties[2];
						double length_ra = properties[3];

						// Increase number of branches
						if(length != 0)
						{
							if(this.auxPoint == null)
								this.auxPoint = nextPoint;

							this.numberOfBranches[iTree]++;

							// Initial vertex
							Vertex initialVertex = null;
							for(int k = 0; k < this.junctionVertex[iTree].length; k++)
								if(this.junctionVertex[iTree][k].isVertexPoint(junctionCoord))
								{
									initialVertex = this.junctionVertex[iTree][k];
									break;
								}


							// If the final point is a slab, then we add the path to the
							// neighbor junction voxel not belonging to the initial vertex
							// (unless it is a self loop)
							if(isSlab(this.auxPoint))
							{
								final Point aux = this.auxPoint;
								//IJ.log("Looking for " + this.auxPoint + " in the list of vertices...");
								this.auxPoint = getVisitedJunctionNeighbor(this.auxPoint, initialVertex);
								this.auxFinalVertex = findPointVertex(this.junctionVertex[iTree], this.auxPoint);
								if(this.auxPoint == null)
								{
									//IJ.log("Point "+ aux + " has not neighbor end junction! (inner loop)");
									// Inner loop
									this.auxFinalVertex = initialVertex;
									this.auxPoint = aux;
								}
								length += calculateDistance(this.auxPoint, aux);
							}

							if(debug)
								IJ.log("increased number of branches, length = " + length + " (last point = " + this.auxPoint + ")");
							// update maximum branch length
							if(length > this.maximumBranchLength[iTree])
							{
								this.maximumBranchLength[iTree] = length;
							}

							// Add the distance between the main vertex of the junction 
							// and the initial junction vertex of the branch (this prevents from
							// having branches in the graph larger than the calculated branch length)
							length += calculateDistance(initialVertex.getPoints().get(0), junctionCoord);

							// Create graph branch
							// Add branch to graph
							if(debug)
								IJ.log("adding branch from " + initialVertex.getPoints().get(0) + " to " + this.auxFinalVertex.getPoints().get(0));
							this.graph[iTree].addEdge(new Edge(initialVertex, this.auxFinalVertex, this.slabList, length, color3rd, color, length_ra));
							// Increase total length of branches
							branchLength += length;

							// update maximum branch length
							if(length > this.maximumBranchLength[iTree])
							{
								this.maximumBranchLength[iTree] = length;
							}
						}
					}
					else
						setVisited(nextPoint, true);

					nextPoint = getNextUnvisitedVoxel(junctionCoord);
				}
			}
		}

		if(debug)
			IJ.log( " --------------------------- ");

		// Finally visit branches starting at slabs (special case for circular trees)
		if(this.startingSlabTree[iTree].size() == 1)
		{
			if(debug)
				IJ.log("visit from slabs");

			final Point startCoord = this.startingSlabTree[iTree].get(0);

			// Create circular graph (only one vertex)
			final Vertex v1 = new Vertex();
			v1.addPoint(startCoord);
			this.graph[iTree].addVertex(v1);


			this.slabList = new ArrayList<Point>();
			this.slabList.add(startCoord);

			this.numberOfSlabs[iTree]++;

			// visit branch until finding visited voxel.
			double[] properties = visitBranch(startCoord, iTree);
			final double length = properties[0];
			double color3rd = properties[1];
			double color = properties[2];
			double length_ra = properties[3];

			if(length != 0)
			{
				// increase number of branches
				this.numberOfBranches[iTree]++;
				branchLength += length;

				// update maximum branch length
				if(length > this.maximumBranchLength[iTree])
				{
					this.maximumBranchLength[iTree] = length;
				}
			}

			// Create circular edge
			this.graph[iTree].addEdge(new Edge(v1, v1, this.slabList, length, color3rd, color, length_ra));
		}

		if(debug)
			IJ.log( " --------------------------- ");

		if(this.numberOfBranches[iTree] == 0)
			return;
		// Average length
		this.averageBranchLength[iTree] = branchLength / this.numberOfBranches[iTree];

		if(debug)
		{
			IJ.log("Num of vertices = " + this.graph[iTree].getVertices().size() + " num of edges = " + this.graph[iTree].getEdges().size());
			for(int i = 0; i < this.graph[iTree].getVertices().size(); i++)
			{
				Vertex v = this.graph[iTree].getVertices().get(i);
				IJ.log(" vertex " + v.getPoints().get(0) + " has neighbors: ");
				for(int j = 0; j < v.getBranches().size(); j++)
				{
					final Vertex v1 = v.getBranches().get(j).getV1();
					final Vertex oppositeVertex = v1.equals(v) ? v.getBranches().get(j).getV2() : v1 ;
					IJ.log(j + ": " + oppositeVertex.getPoints().get(0));
				}

			}

			IJ.log( " --------------------------- ");
			for(int i = 0; i < this.junctionVertex[iTree].length; i ++)
			{
				IJ.log("Junction #" + i + " is formed by: ");
				for(int j = 0; j < this.junctionVertex[iTree][i].getPoints().size(); j++)
					IJ.log(j + ": " + this.junctionVertex[iTree][i].getPoints().get(j));
			}
		}

	} // end visitSkeleton

	/* -----------------------------------------------------------------------*/
	/**
	 * Color the different trees in the skeleton.
	 * 
	 * @param taggedImage
	 * 
	 * @return image with every tree tagged with a different number 
	 */
	private ImageStack markTrees(ImageStack taggedImage) 
	{
		if(debug)
			IJ.log("=== Mark Trees ===");
		// Create output image
		ImageStack outputImage = new ImageStack( this.width, this.height );
		for (int z = 0; z < depth; z++)
		{
			outputImage.addSlice(taggedImage.getSliceLabel(z+1), new FloatProcessor(this.width, this.height));
		}

		this.numOfTrees = 0;

		int color = 0;

		// Visit trees starting at end points
		for(int i = 0; i < this.totalNumberOfEndPoints; i++)
		{
			Point endPointCoord = this.listOfEndPoints.get(i);

			if(isVisited(endPointCoord))
				continue;

			color++;

			if(color == Integer.MAX_VALUE)
			{
				IJ.error("More than " + (Integer.MAX_VALUE-1) +
						" skeletons in the image. AnalyzeSkeleton can only process up to "+ (Integer.MAX_VALUE-1));
				return null;
			}

			if(debug)
				IJ.log("-- Visit tree from end-point:");
			// Visit the entire tree.
			int numOfVoxelsInTree = visitTree(endPointCoord, outputImage, color);

			// increase number of trees
			this.numOfTrees++;
		}

		// Visit trees starting at junction points 
		// (some circular trees do not have end points)
		// Visit trees starting at end points
		for(int i = 0; i < this.totalNumberOfJunctionVoxels; i++)
		{
			Point junctionCoord = this.listOfJunctionVoxels.get(i);
			if(isVisited(junctionCoord))
				continue;

			color++;

			if(color == Integer.MAX_VALUE)
			{
				IJ.error("More than " + (Integer.MAX_VALUE-1) +
						" skeletons in the image. AnalyzeSkeleton can only process up to "+ (Integer.MAX_VALUE-1));
				return null;
			}

			if(debug)
				IJ.log("-- Visit tree from junction:");

			// else, visit branch until next junction or end point.
			int length = visitTree(junctionCoord, outputImage, color);

			if(length == 0)
			{
				color--; // the color was not used
				continue;
			}

			// increase number of trees
			this.numOfTrees++;
		}

		// Check for unvisited slab voxels
		// (just in case there are circular trees without junctions)
		for(int i = 0; i < this.listOfSlabVoxels.size(); i++)
		{
			Point p = (Point) this.listOfSlabVoxels.get(i);
			if(isVisited(p) == false)
			{
				// Mark that voxel as the start point of the circular skeleton
				this.listOfStartingSlabVoxels.add(p);

				color++;

				if(color == Integer.MAX_VALUE)
				{
					IJ.error("More than " + (Integer.MAX_VALUE-1) +
							" skeletons in the image. AnalyzeSkeleton can only process up to "+ (Integer.MAX_VALUE-1));
					return null;
				}

				if(debug)
					IJ.log("-- Visit tree from slab:");

				// else, visit branch until next junction or end point.
				int length = visitTree(p, outputImage, color);

				if(length == 0)
				{
					color--; // the color was not used
					continue;
				}

				// increase number of trees
				this.numOfTrees++;
			}
		}



		//System.out.println("Number of trees = " + this.numOfTrees);

		// Show tree image.
		if(debug)
		{
			ImagePlus treesIP = new ImagePlus("Trees skeleton", outputImage);
			treesIP.show();

			// Set same calibration as the input image
			treesIP.setCalibration(this.imRef.getCalibration());

			// We apply the Fire LUT and reset the min and max to be between 0-255.
			IJ.run("Fire");

			//IJ.resetMinAndMax();
			treesIP.resetDisplayRange();
			treesIP.updateAndDraw();
		}

		// Reset visited variable
		resetVisited();

		//IJ.log("Number of trees: " + this.numOfTrees + ", # colors = " + color);

		return outputImage;

	} /* end markTrees */


	// --------------------------------------------------------------
	/**
	 * Visit tree marking the voxels with a reference tree color.
	 * 
	 * @param startingPoint starting tree point
	 * @param outputImage 3D image to visit
	 * @param color reference tree color
	 * @return number of voxels in the tree
	 */
	private int visitTree(Point startingPoint, ImageStack outputImage,
			int color) 
	{
		int numOfVoxels = 0;

		if(debug)
			IJ.log("visiting " + startingPoint + " color = " + color);

		if(isVisited(startingPoint))
			return 0;

		// Set pixel color
		this.setPixel(outputImage, startingPoint.x, startingPoint.y, startingPoint.z, color);
		setVisited(startingPoint, true);

		ArrayList <Point> toRevisit = new ArrayList <Point>();

		// Add starting point to revisit list if it is a junction
		if(isJunction(startingPoint))
			toRevisit.add(startingPoint);

		Point nextPoint = getNextUnvisitedVoxel(startingPoint);

		while(nextPoint != null || toRevisit.size() != 0)
		{
			if(nextPoint != null)
			{
				if(!isVisited(nextPoint))
				{
					numOfVoxels++;
					if(debug)
						IJ.log("visiting " + nextPoint+ " color = " + color);

					// Set color and visit flat
					this.setPixel(outputImage, nextPoint.x, nextPoint.y, nextPoint.z, color);
					setVisited(nextPoint, true);

					// If it is a junction, add it to the revisit list
					if(isJunction(nextPoint))
						toRevisit.add(nextPoint);

					// Calculate next point to visit
					nextPoint = getNextUnvisitedVoxel(nextPoint);
				}
			}
			else // revisit list
			{
				nextPoint = toRevisit.get(0);
				if(debug)
					IJ.log("revisiting " + nextPoint+ " color = " + color);

				// Calculate next point to visit
				nextPoint = getNextUnvisitedVoxel(nextPoint);
				// Maintain junction in the list until there is no more branches
				if (nextPoint == null)
					toRevisit.remove(0);
			}
		}

		return numOfVoxels;
	} // end method visitTree

	// -----------------------------------------------------------------------
	/**
	 * Visit a branch and calculate length.
	 * 
	 * @param startingPoint starting coordinates
	 * @return branch length
	 * 
	 * @deprecated
	 */
	private double visitBranch(Point startingPoint) 
	{
		double length = 0;

		// mark starting point as visited
		setVisited(startingPoint, true);

		// Get next unvisited voxel
		Point nextPoint = getNextUnvisitedVoxel(startingPoint);

		if (nextPoint == null)
			return 0;

		Point previousPoint = startingPoint;

		// We visit the branch until we find an end point or a junction
		while(nextPoint != null && isSlab(nextPoint))
		{
			// Add length
			length += calculateDistance(previousPoint, nextPoint);

			// Mark as visited
			setVisited(nextPoint, true);

			// Move in the graph
			previousPoint = nextPoint;
			nextPoint = getNextUnvisitedVoxel(previousPoint);
		}


		if(nextPoint != null)
		{
			// Add distance to last point
			length += calculateDistance(previousPoint, nextPoint);

			// Mark last point as visited
			setVisited(nextPoint, true);
		}

		this.auxPoint = previousPoint;

		return length;
	}// end visitBranch

	// -----------------------------------------------------------------------
	/**
	 * Visit a branch and calculate length in a specific tree
	 * 
	 * @param startingPoint starting coordinates
	 * @param iTree tree index
	 * @return branch length and color
	 */
	private double[] visitBranch(Point startingPoint, int iTree) 
	{
		//IJ.log("startingPoint = (" + startingPoint.x + ", " + startingPoint.y + ", " + startingPoint.z + ")");
		double length = 0;
		double intensity = 0.0;
		double intensity3rd = 0.0;
		double length_ra = 0.0;
		double[] ret = new double[4];

		List<Point> pointHistory= new ArrayList<Point>(0);
		pointHistory.add(0,startingPoint);


		// mark starting point as visited
		setVisited(startingPoint, true);

		// Get next unvisited voxel
		Point nextPoint = getNextUnvisitedVoxel(startingPoint);

		if (nextPoint == null)
			return ret;

		Point previousPoint = startingPoint;

		// We visit the branch until we find an end point or a junction
		while(nextPoint != null && isSlab(nextPoint))
		{
			this.numberOfSlabs[iTree]++;

			// Add slab voxel to the edge
			this.slabList.add(nextPoint);

			// Add length
			length += calculateDistance(previousPoint, nextPoint);
			pointHistory.add(0,nextPoint);

			length_ra += calculateDistance(pointHistory);

			// Mark as visited
			setVisited(nextPoint, true);

			// Move in the graph
			previousPoint = nextPoint;
			if (debug)
				IJ.log("visiting " + previousPoint);
			nextPoint = getNextUnvisitedVoxel(previousPoint);
		}

		// If we find an unvisited end-point or junction, we set it
		// as final vertex of the branch
		if(nextPoint != null)
		{
			// Add distance to last point
			length += calculateDistance(previousPoint, nextPoint);
			pointHistory.add(0,nextPoint);
			length_ra += calculateDistance(pointHistory);


			// Mark last point as visited
			setVisited(nextPoint, true);

			// Mark final vertex
			if(isEndPoint(nextPoint))
			{
				if(debug)
					IJ.log("found unvisited end point: " + nextPoint);
				this.auxFinalVertex = new Vertex();
				this.auxFinalVertex.addPoint(nextPoint);
			}
			else if(isJunction(nextPoint))
			{
				if(debug)
					IJ.log("found unvisited junction point: " + nextPoint);
				this.auxFinalVertex = findPointVertex(this.junctionVertex[iTree], nextPoint);
				// Add the length to the first point of the vertex (to prevent later from having
				// euclidean distances larger than the actual distance)
				length += calculateDistance(auxFinalVertex.getPoints().get(0), nextPoint);
				length_ra += calculateDistance(auxFinalVertex.getPoints().get(0), nextPoint);
				/*
				int j = 0;
				for(j = 0; j < this.junctionVertex[iTree].length; j++)
					if(this.junctionVertex[iTree][j].isVertexPoint(nextPoint))
					{
						this.auxFinalVertex = this.junctionVertex[iTree][j];
						IJ.log(" " + nextPoint + " belongs to junction " + this.auxFinalVertex.getPoints().get(0));
						break;
					}
				if(j == this.junctionVertex[iTree].length)
					IJ.log("point " + nextPoint + " was not found in vertex list!");
					*/
			}

			this.auxPoint = nextPoint;
		}
		else
			this.auxPoint = previousPoint;

		//IJ.log("finalPoint = (" + nextPoint.x + ", " + nextPoint.y + ", " + nextPoint.z + ")");

		// calculate average intensity (thickness) value, but only take the inner third of a branch.
		// at both ends the intensity (thickness) is most likely affected by junctions.


		int size = pointHistory.size();
		int start = (int)(size/3.0);
		int end = (int)(2*size/3.0);

		for (int i = 0; i < size ;i++){
			int value = getPixel(this.inputImage, pointHistory.get(i));
			if (value < 0){
				value += 256;
			}
			intensity += value;
			if (i>=start && i<end){
				intensity3rd += value;
			}
		}

		intensity /= size; 
		intensity3rd /= (end-start);

		ret[0] = length;
		ret[1] = intensity3rd;
		ret[2] = intensity;
		ret[3] = length_ra;

		return ret;
	} // end visitBranch

	// -----------------------------------------------------------------------
	/**
	 * Find vertex in an array given a specific vertex point.
	 * 
	 * @param vertex array of search
	 * @param p vertex point
	 * @return vertex containing that point
	 */
	public Vertex findPointVertex(Vertex[] vertex, Point p)
	{
		int j = 0;
		for(j = 0; j < vertex.length; j++)
			if(vertex[j].isVertexPoint(p))
			{
				if(debug)
					IJ.log(" " + p + " belongs to junction " + vertex[j].getPoints().get(0));
				return vertex[j];
			}
		if(debug)
			IJ.log("point " + p + " was not found in vertex list! (vertex.length= " + vertex.length +")");
		return null;
	}

	// -----------------------------------------------------------------------
	/**
	 * Calculate Euclidean distance between two points in 3D.
	 * 
	 * @param point1 first point coordinates
	 * @param point2 second point coordinates
	 * @return distance (in the corresponding units)
	 */
	public double calculateDistance(Point point1, Point point2) 
	{
		final double dx = (point1.x - point2.x) * this.imRef.getCalibration().pixelWidth;
		final double dy = (point1.y - point2.y) * this.imRef.getCalibration().pixelHeight;
		final double dz = (point1.z - point2.z) * this.imRef.getCalibration().pixelDepth;
		return Math.sqrt(dx * dx + dy * dy + dz * dz);
	}

	// -----------------------------------------------------------------------
		/**
		 * Calculate linear corrected distance between two points in 3D.
		 * Uses the PointHistory of a branch and its 5 last points.
		 * 
		 * @param Points - the last visited Points (most recent has Index 0)
		 * @return linear corrected distance between the last two Points (in the corresponding units)
		 */
	private double calculateDistance(List<Point> Points)
	{
		int indexOfLast = Points.size()-1;

		//no Distance to be calculated here...
		if (indexOfLast < 1) return 0;

		// Point Of Interest
		int poi = 5;

		//poi is indexOflast if List is shorter than 5
		if (indexOfLast < 5){
			poi = indexOfLast;
		}
		return calculateDistance(Points.get(poi), Points.get(0))/poi;
	}

	// -----------------------------------------------------------------------
	/**
	 * Calculate number of junction skipping neighbor junction voxels
	 * 
	 * @param treeIS tree stack
	 */
	private void groupJunctions(ImageStack treeIS) 
	{
		// Mark all unvisited
		resetVisited();

		for (int iTree = 0; iTree < this.numOfTrees; iTree++)
		{
			// Visit list of junction voxels
			for(int i = 0; i < this.numberOfJunctionVoxels[iTree]; i ++)
			{
				Point pi = this.junctionVoxelTree[iTree].get(i);

				if(! isVisited(pi))
					fusionNeighborJunction(pi, this.listOfSingleJunctions[iTree]);
			}
		}

		// Count number of single junctions for every tree in the image
		for (int iTree = 0; iTree < this.numOfTrees; iTree++)
		{
			if(debug)
				IJ.log("this.listOfSingleJunctions["+iTree+"].size() = " + this.listOfSingleJunctions[iTree].size());

			this.numberOfJunctions[iTree] = this.listOfSingleJunctions[iTree].size();

			// Create array of junction vertices for the graph
			this.junctionVertex[iTree] = new Vertex[this.listOfSingleJunctions[iTree].size()];

			for(int j = 0 ; j < this.listOfSingleJunctions[iTree].size(); j++)
			{
				final ArrayList<Point> list = this.listOfSingleJunctions[iTree].get(j);
				this.junctionVertex[iTree][j] = new Vertex();
				for(final Point p : list)
					this.junctionVertex[iTree][j].addPoint(p);

			}
		}

		// Mark all unvisited
		resetVisited();
	}

	// -----------------------------------------------------------------------
	/**
	 * Reset visit variable and set it to false.
	 */
	private void resetVisited()
	{
		// Reset visited variable
		this.visited = null;
		this.visited = new boolean[this.width][this.height][this.depth];
		/*
		for(int i = 0; i < this.width; i ++)
			for(int j = 0; j < this.height; j++)
				for(int k = 0; k < this.depth; k++)
					this.visited[i][j][k] = false;
		*/
	}

	// -----------------------------------------------------------------------
	/**
	 * Fusion neighbor junctions voxels into the same list.
	 * 
	 * @param startingPoint starting junction voxel
	 * @param singleJunctionsList list of single junctions
	 */
	private void fusionNeighborJunction(Point startingPoint,
			ArrayList<ArrayList<Point>> singleJunctionsList) 
	{
		// Create new group of junctions
		ArrayList <Point> newGroup = new ArrayList<Point>();
		newGroup.add(startingPoint);

		// Mark the starting junction as visited
		setVisited(startingPoint, true);

		// Look for neighbor junctions and add them to the new group
		ArrayList <Point> toRevisit = new ArrayList <Point>();
		toRevisit.add(startingPoint);

		Point nextPoint = getNextUnvisitedJunctionVoxel(startingPoint);

		while(nextPoint != null || toRevisit.size() != 0)
		{
			if(nextPoint != null && !isVisited(nextPoint))
			{
				// Add to the group
				newGroup.add(nextPoint);
				// Mark as visited
				setVisited(nextPoint, true);

				// add it to the revisit list
				toRevisit.add(nextPoint);

				// Calculate next junction point to visit
				nextPoint = getNextUnvisitedJunctionVoxel(nextPoint);
			}
			else // revisit list
			{
				nextPoint = toRevisit.get(0);
				//IJ.log("visiting " + nextPoint + " color = " + color);

				// Calculate next point to visit
				nextPoint = getNextUnvisitedJunctionVoxel(nextPoint);
				// Maintain junction in the list until there is no more branches
				if (nextPoint == null)
					toRevisit.remove(0);
			}
		}

		// Add group to the single junction list
		singleJunctionsList.add(newGroup);

	}// end method fusionNeighborJunction

	// -----------------------------------------------------------------------
	/**
	 *  Check if two groups of voxels are neighbors.
	 *  
	 * @param g1 first group
	 * @param g2 second group
	 * 
	 * @return true if the groups have any neighbor voxel
	 */
	boolean checkNeighborGroups(ArrayList <Point> g1,  ArrayList <Point> g2)
	{
		for(int i =0 ; i  < g1.size(); i++)
		{
			Point pi = g1.get(i);
			for(int j =0 ; j  < g2.size(); j++)
			{
				Point pj = g2.get(j);
				if(isNeighbor(pi, pj))
					return true;
			}
		}
		return false;
	}


	// -----------------------------------------------------------------------
	/**
	 * Calculate number of triple and quadruple points in the skeleton. Triple and 
	 * quadruple points are junctions with exactly 3 and 4 branches respectively.
	 */
	private void calculateTripleAndQuadruplePoints() 
	{
		for (int iTree = 0; iTree < this.numOfTrees; iTree++)
		{
			// Visit the groups of junction voxels
			for(int i = 0; i < this.numberOfJunctions[iTree]; i ++)
			{

				ArrayList <Point> groupOfJunctions = this.listOfSingleJunctions[iTree].get(i);

				// Count the number of slab and end-points neighbors of every voxel in the group
				int nBranch = 0;
				for(int j = 0; j < groupOfJunctions.size(); j++)
				{
					Point pj = groupOfJunctions.get(j);

					// Get neighbors and check the slabs or end-points
					byte[] neighborhood = this.getNeighborhood(this.taggedImage, pj.x, pj.y, pj.z);
					for(int k = 0; k < 27; k++)
						if (neighborhood[k] == AnalyzeSkeleton_.SLAB 
							|| neighborhood[k] == AnalyzeSkeleton_.END_POINT)
							nBranch++;
				}
				// If the junction has only 3 slab/end-point neighbors, then it is a triple point
				if (nBranch == 3)
					this.numberOfTriplePoints[iTree] ++;
				else if(nBranch == 4) // quadruple point if 4
					this.numberOfQuadruplePoints[iTree] ++;
			}

		}

	}// end calculateTripleAndQuadruplePoints


	/* -----------------------------------------------------------------------*/
	/**
	 * Calculate if two points are neighbors.
	 * 
	 * @param point1 first point
	 * @param point2 second point
	 * @return true if the points are neighbors (26-pixel neighborhood)
	 */
	private boolean isNeighbor(Point point1, Point point2) 
	{
		return Math.sqrt(  Math.pow( (point1.x - point2.x), 2) 
		          + Math.pow( (point1.y - point2.y), 2)
		          + Math.pow( (point1.z - point2.z), 2)) <= Math.sqrt(3);
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if the point is slab.
	 *  
	 * @param point actual point
	 * @return true if the point has slab status
	 */
	private boolean isSlab(Point point) 
	{
		return getPixel(this.taggedImage, point.x, point.y, point.z) == AnalyzeSkeleton_.SLAB;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if the point is a junction.
	 *  
	 * @param point actual point
	 * @return true if the point has slab status
	 */
	private boolean isJunction(Point point) 
	{
		return getPixel(this.taggedImage, point.x, point.y, point.z) == AnalyzeSkeleton_.JUNCTION;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if the point is an end point.
	 *  
	 * @param point actual point
	 * @return true if the point has slab status
	 */
	private boolean isEndPoint(Point point) 
	{
		return getPixel(this.taggedImage, point.x, point.y, point.z) == AnalyzeSkeleton_.END_POINT;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if the point is an 'unprotected' end point.
	 *
	 * @param point actual point
	 * @param roi Only points outside this ROI will have end-point status. ROI
	 *            can be associated to a single image in the stack (or all) as
	 *            per {@link ij.gui.Roi#getPosition ij.gui.Roi.getPosition()}
	 * @return true if the point has end-point status
	 */
	private boolean isEndPoint(Point point, Roi roi)
	{
		boolean endPoint = isEndPoint(point);
		if (endPoint && roi != null)
		{
			// Is roi associated with all the images in the ImageStack?
			boolean roiContainsZ = (roi.getPosition()==0) ? true : roi.getPosition()==point.z;

			// Maintain end-point status only if point lies outside roi boundaries
			endPoint = !(roi.contains(point.x, point.y) && roiContainsZ);
		}
		return endPoint;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if the point is a junction.
	 *  
	 * @param x x- voxel coordinate
	 * @param y y- voxel coordinate
	 * @param z z- voxel coordinate
	 * @return true if the point has slab status
	 */
	private boolean isJunction(int x, int y, int z) 
	{
		return getPixel(this.taggedImage, x, y, z) == AnalyzeSkeleton_.JUNCTION;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Get next unvisited neighbor voxel.
	 * 
	 * @param point starting point
	 * @return unvisited neighbor or null if all neighbors are visited
	 */
	private Point getNextUnvisitedVoxel(Point point) 
	{
		Point unvisitedNeighbor = null;

		// Check neighbors status
		for(int x = -1; x < 2; x++)
			for(int y = -1; y < 2; y++)
				for(int z = -1; z < 2; z++)
				{
					if(x == 0 && y == 0 && z == 0)
						continue;

					if(getPixel(this.inputImage, point.x + x, point.y + y, point.z + z) != 0
						&& isVisited(point.x + x, point.y + y, point.z + z) == false)
					{
						unvisitedNeighbor = new Point(point.x + x, point.y + y, point.z + z);
						break;
					}

				}

		return unvisitedNeighbor;
	}// end getNextUnvisitedVoxel

	/* -----------------------------------------------------------------------*/
	/**
	 * Get next unvisited junction neighbor voxel.
	 * 
	 * @param point starting point
	 * @return unvisited neighbor or null if all neighbors are visited
	 */
	private Point getNextUnvisitedJunctionVoxel(Point point) 
	{
		Point unvisitedNeighbor = null;

		// Check neighbors status
		for(int x = -1; x < 2; x++)
			for(int y = -1; y < 2; y++)
				for(int z = -1; z < 2; z++)
				{
					if(x == 0 && y == 0 && z == 0)
						continue;

					if(getPixel(this.inputImage, point.x + x, point.y + y, point.z + z) != 0
						&& isVisited(point.x + x, point.y + y, point.z + z) == false 
						&& isJunction(point.x + x, point.y + y, point.z + z))
					{
						unvisitedNeighbor = new Point(point.x + x, point.y + y, point.z + z);
						break;
					}

				}

		return unvisitedNeighbor;
	}// end getNextUnvisitedJunctionVoxel 

	// -----------------------------------------------------------------------
	/**
	 * Get next visited junction neighbor voxel excluding the ones belonging
	 * to a give vertex
	 * 
	 * @param point starting point
	 * @param exclude exclusion vertex
	 * @return unvisited neighbor or null if all neighbors are visited
	 */
	private Point getVisitedJunctionNeighbor(Point point, Vertex exclude) 
	{
		Point finalNeighbor = null;

		// Check neighbors status
		for(int x = -1; x < 2; x++)
			for(int y = -1; y < 2; y++)
				for(int z = -1; z < 2; z++)
				{
					if(x == 0 && y == 0 && z == 0)
						continue;

					final Point neighbor = new Point( point.x + x, point.y + y, point.z + z);

					if(getPixel(this.inputImage, neighbor) != 0
						&& isVisited(neighbor)
						&& isJunction(neighbor)
						&& ! exclude.getPoints().contains(neighbor))
					{
						finalNeighbor = neighbor;
						break;
					}

				}

		return finalNeighbor;
	}// end getNextUnvisitedJunctionVoxel 

	// -----------------------------------------------------------------------
	/**
	 * Check if a voxel is visited taking into account the borders. 
	 * Out of range voxels are considered as visited. 
	 * 
	 * @param point
	 * @return true if the voxel is visited
	 */
	private boolean isVisited(Point point) 
	{
		return isVisited(point.x, point.y, point.z);
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Check if a voxel is visited taking into account the borders. 
	 * Out of range voxels are considered as visited. 
	 * 
	 * @param x x- voxel coordinate
	 * @param y y- voxel coordinate
	 * @param z z- voxel coordinate
	 * @return true if the voxel is visited
	 */
	private boolean isVisited(int x, int y, int z) 
	{
		if(x >= 0 && x < this.width && y >= 0 && y < this.height && z >= 0 && z < this.depth)
			return this.visited[x][y][z];
		return true;
	}


	/* -----------------------------------------------------------------------*/
	/**
	 * Set value in the visited flags matrix.
	 * 
	 * @param x x- voxel coordinate
	 * @param y y- voxel coordinate
	 * @param z z- voxel coordinate
	 * @param b
	 */
	private void setVisited(int x, int y, int z, boolean b) 
	{
		if(x >= 0 && x < this.width && y >= 0 && y < this.height && z >= 0 && z < this.depth)
			this.visited[x][y][z] = b;
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Set value in the visited flags matrix.
	 * 
	 * @param point voxel coordinates
	 * @param b visited flag value
	 */
	private void setVisited(Point point, boolean b) 
	{
		int x = point.x;
		int y = point.y;
		int z = point.z;

 		setVisited(x, y, z, b);
	}

	/* -----------------------------------------------------------------------*/
	/**
	 * Tag skeleton dividing the voxels between end points, junctions and slabs.
	 *  
	 * @param inputImage2 skeleton image to be tagged
	 * @return tagged skeleton image
	 */
	private ImageStack tagImage(ImageStack inputImage2) 
	{
		// Create output image
		ImageStack outputImage = new ImageStack(this.width, this.height, inputImage2.getColorModel());

		// Tag voxels
		for (int z = 0; z < depth; z++)
		{
			outputImage.addSlice(inputImage2.getSliceLabel(z+1), new ByteProcessor(this.width, this.height));
			for (int x = 0; x < width; x++) 
				for (int y = 0; y < height; y++)
				{
					if(getPixel(inputImage2, x, y, z) != 0)
					{
						int numOfNeighbors = getNumberOfNeighbors(inputImage2, x, y, z);
						if(numOfNeighbors < 2)
						{
							setPixel(outputImage, x, y, z, AnalyzeSkeleton_.END_POINT);
							this.totalNumberOfEndPoints++;
							Point endPoint = new Point(x, y, z);
							this.listOfEndPoints.add(endPoint);
						}
						else if(numOfNeighbors > 2)
						{
							setPixel(outputImage, x, y, z, AnalyzeSkeleton_.JUNCTION);
							Point junction = new Point(x, y, z);
							this.listOfJunctionVoxels.add(junction);
							this.totalNumberOfJunctionVoxels++;
						}
						else
						{
							setPixel(outputImage, x, y, z, AnalyzeSkeleton_.SLAB);
							Point slab = new Point(x, y, z);
							this.listOfSlabVoxels.add(slab);
							this.totalNumberOfSlabs++;
						}
					}
				}
		}

		return outputImage;
	}// end method tagImage 

	/* -----------------------------------------------------------------------*/
	/**
	 * Get number of neighbors of a voxel in a 3D image (0 border conditions).
	 * 
	 * @param image 3D image (ImageStack)
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding 27-pixels neighborhood (0 if out of image)
	 */
	private int getNumberOfNeighbors(ImageStack image, int x, int y, int z)
	{
		int n = 0;
		byte[] neighborhood = getNeighborhood(image, x, y, z);

		for(int i = 0; i < 27; i ++)
			if(neighborhood[i] != 0)
				n++;
		// We return n-1 because neighborhood includes the actual voxel.
		return (n-1);
	}// end method getNumberOfNeighbors

	// -----------------------------------------------------------------------
	/**
	 * Get average 3x3x3 neighborhood pixel value of a given point
	 * 
	 * @param image input image
	 * @param p image coordinates
	 */
	private double getAverageNeighborhoodValue(ImageStack image, Point p)
	{
		byte[] neighborhood = getNeighborhood(image, p);

		double avg = 0;
		for(int i = 0; i < neighborhood.length; i++)
			avg += (double) (neighborhood[i] & 0xFF);
		if(neighborhood.length > 0)
			return avg / (double) neighborhood.length;
		else
			return 0;
	}// end method getAverageNeighborhoodValue

	// -----------------------------------------------------------------------
	/**
	 * Get average neighborhood pixel value of a given point.
	 * 
	 * @param image input image
	 * @param p image coordinates
	 * @param x_offset x- neighborhood offset
	 * @param y_offset y- neighborhood offset
	 * @param z_offset z- neighborhood offset
	 * @return average neighborhood pixel value 
	 */
	public static double getAverageNeighborhoodValue(
			final ImageStack image, 
			final Point p,
			final int x_offset,
			final int y_offset,
			final int z_offset)
	{
		byte[] neighborhood = getNeighborhood(image, p, x_offset, y_offset, z_offset);

		double avg = 0;
		for(int i = 0; i < neighborhood.length; i++)
			avg += (double) (neighborhood[i] & 0xFF);
		if(neighborhood.length > 0)
			return avg / (double) neighborhood.length;
		else
			return 0;
	}// end method getAverageNeighborhoodValue


	// -----------------------------------------------------------------------
	/**
	 * Get neighborhood of a pixel in a 3D image (0 border conditions).
	 * 
	 * @param image 3D image (ImageStack)
	 * @param p point coordinates
	 * @param x_offset x- neighborhood offset
	 * @param y_offset y- neighborhood offset
	 * @param z_offset z- neighborhood offset
	 * @return corresponding neighborhood (0 if out of image)
	 */
	public static byte[] getNeighborhood(
			final ImageStack image, 
			final Point p, 
			final int x_offset, 
			final int y_offset, 
			final int z_offset)
	{
		final byte[] neighborhood = new byte[(2*x_offset+1) * (2*y_offset+1) * (2*z_offset+1)];

		for(int l= 0, k = p.z - z_offset; k <= p.z + z_offset; k++)
			for(int j = p.y - y_offset; j <= p.y + y_offset; j++)
				for(int i = p.x - x_offset; i <= p.x + x_offset; i++, l++)
					neighborhood[l] = getPixel(image, i,   j,   k);
		return neighborhood;
	} // end getNeighborhood 

	// -----------------------------------------------------------------------
	/**
	 * Get neighborhood of a pixel in a 3D image (0 border conditions).
	 * 
	 * @param image 3D image (ImageStack)
	 * @param p 3D point coordinates
	 * @return corresponding 27-pixels neighborhood (0 if out of image)
	 */
	private byte[] getNeighborhood(ImageStack image, Point p)
	{
		return getNeighborhood(image, p.x, p.y, p.z);
	}

	// -----------------------------------------------------------------------
	/**
	 * Get neighborhood of a pixel in a 3D image (0 border conditions).
	 * 
	 * @param image 3D image (ImageStack)
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding 27-pixels neighborhood (0 if out of image)
	 */
	private byte[] getNeighborhood(ImageStack image, int x, int y, int z)
	{
		byte[] neighborhood = new byte[27];

		neighborhood[ 0] = getPixel(image, x-1, y-1, z-1);
		neighborhood[ 1] = getPixel(image, x  , y-1, z-1);
		neighborhood[ 2] = getPixel(image, x+1, y-1, z-1);

		neighborhood[ 3] = getPixel(image, x-1, y,   z-1);
		neighborhood[ 4] = getPixel(image, x,   y,   z-1);
		neighborhood[ 5] = getPixel(image, x+1, y,   z-1);

		neighborhood[ 6] = getPixel(image, x-1, y+1, z-1);
		neighborhood[ 7] = getPixel(image, x,   y+1, z-1);
		neighborhood[ 8] = getPixel(image, x+1, y+1, z-1);

		neighborhood[ 9] = getPixel(image, x-1, y-1, z  );
		neighborhood[10] = getPixel(image, x,   y-1, z  );
		neighborhood[11] = getPixel(image, x+1, y-1, z  );

		neighborhood[12] = getPixel(image, x-1, y,   z  );
		neighborhood[13] = getPixel(image, x,   y,   z  );
		neighborhood[14] = getPixel(image, x+1, y,   z  );

		neighborhood[15] = getPixel(image, x-1, y+1, z  );
		neighborhood[16] = getPixel(image, x,   y+1, z  );
		neighborhood[17] = getPixel(image, x+1, y+1, z  );

		neighborhood[18] = getPixel(image, x-1, y-1, z+1);
		neighborhood[19] = getPixel(image, x,   y-1, z+1);
		neighborhood[20] = getPixel(image, x+1, y-1, z+1);

		neighborhood[21] = getPixel(image, x-1, y,   z+1);
		neighborhood[22] = getPixel(image, x,   y,   z+1);
		neighborhood[23] = getPixel(image, x+1, y,   z+1);

		neighborhood[24] = getPixel(image, x-1, y+1, z+1);
		neighborhood[25] = getPixel(image, x,   y+1, z+1);
		neighborhood[26] = getPixel(image, x+1, y+1, z+1);

		return neighborhood;
	} // end getNeighborhood 


	// -----------------------------------------------------------------------
	/**
	 * Get pixel in 3D image (0 border conditions) 
	 * 
	 * @param image 3D image
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding pixel (0 if out of image)
	 */
	public static byte getPixel(final ImageStack image, final int x, final int y, final int z)
	{
		final int width = image.getWidth();
		final int height = image.getHeight();
		final int depth = image.getSize();

		if(x >= 0 && x < width && y >= 0 && y < height && z >= 0 && z < depth)
			return ((byte[]) image.getPixels(z + 1))[x + y * width];
		else return 0;
	} // end getPixel 

	/* -----------------------------------------------------------------------*/
	/**
	 * Get pixel in 3D image (0 border conditions).
	 * 
	 * @param image 3D image
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @return corresponding pixel (0 if out of image)
	 */
	private float getFloatPixel(
			ImageStack image, 
			int x, 
			int y, 
			int z )
	{
		if(x >= 0 && x < this.width && y >= 0 && y < this.height && z >= 0 && z < this.depth)
			return ( (float[]) image.getPixels(z + 1) )[x + y * this.width];
		else return 0;
	} // end getFloatPixel

	/* -----------------------------------------------------------------------*/
	/**
	 * Get pixel in 3D image (0 border conditions).
	 * 
	 * @param image 3D image
	 * @param point point to be evaluated
	 * @return corresponding pixel (0 if out of image)
	 */
	private float getFloatPixel(
			ImageStack image, 
			Point point )
	{
		return getFloatPixel( image, point.x, point.y, point.z );
	} // end getFloatPixel


	/* -----------------------------------------------------------------------*/
	/**
	 * Get pixel in 3D image (0 border conditions).
	 * 
	 * @param image 3D image
	 * @param point point to be evaluated
	 * @return corresponding pixel (0 if out of image)
	 */
	private byte getPixel(ImageStack image, Point point)
	{
		return getPixel(image, point.x, point.y, point.z);
	} // end getPixel 

	/* -----------------------------------------------------------------------*/
	/**
	 * Set pixel in 3D image.
	 * 
	 * @param image 3D image
	 * @param p point coordinates
	 * @param value pixel value
	 */
	private void setPixel(ImageStack image, Point p, byte value)
	{
		if(p.x >= 0 && p.x < this.width && p.y >= 0 && p.y < this.height && p.z >= 0 && p.z < this.depth)
			((byte[]) image.getPixels(p.z + 1))[p.x + p.y * this.width] = value;
	} // end setPixel 

	/* -----------------------------------------------------------------------*/
	/**
	 * Set pixel in 3D image.
	 * 
	 * @param image 3D image
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @param value pixel value
	 */
	private void setPixel(ImageStack image, int x, int y, int z, byte value)
	{
		if(x >= 0 && x < this.width && y >= 0 && y < this.height && z >= 0 && z < this.depth)
			((byte[]) image.getPixels(z + 1))[x + y * this.width] = value;
	} // end setPixel 

	/* -----------------------------------------------------------------------*/
	/**
	 * Set pixel in 3D (short) image.
	 * 
	 * @param image 3D image
	 * @param x x- coordinate
	 * @param y y- coordinate
	 * @param z z- coordinate (in image stacks the indexes start at 1)
	 * @param value pixel value
	 */
	private void setPixel(ImageStack image, int x, int y, int z, float value)
	{
		if(x >= 0 && x < this.width && y >= 0 && y < this.height && z >= 0 && z < this.depth)
			((float[]) image.getPixels(z + 1))[x + y * this.width] = value;
	} // end setPixel 



	/* -----------------------------------------------------------------------*/
	/**
	 * Show plug-in information.
	 * 
	 */
	void showAbout() 
	{
		IJ.showMessage(
						"About AnalyzeSkeleton...",
						"This plug-in filter analyzes a 2D/3D image skeleton.\n");
	} // end showAbout 


	/**
	 * Determine the longest shortest path using the APSP (all pairs shortest path) 
	 * warshall algorithm
	 * 
	 * @param graph the graph of a tree
	 * @param shortestPathPoints list to store the longest shortest path points
	 * @return longest shortest path length
	 * @author Huub Hovens
	 */
	private double warshallAlgorithm(Graph graph, ArrayList <Point> shortestPathPoints)
	{
		// local fields
		/** vertex 1 of an edge */
		Vertex v1 = null;
		/** vertex 2 of an edge */
		Vertex v2 = null;
		/** the equivalent of row in a matrix */
		int row = 0;
		/** the equivalent of column in a matrix */
		int column = 0;
		/** the value of the longest shortest path */
		double maxPath = 0;
		/** row that contains the longest shortest path value */
		int a = 0;
		/** column that contains the longest shortest path value */
		int b = 0;

		ArrayList< Edge > edgeList = graph.getEdges();
		ArrayList< Vertex > vertexList = graph.getVertices();

		// check for paths of only one vertex
		if( vertexList.size() == 1 )
		{
			shortestPathPoints.add( vertexList.get( 0 ).getPoints().get( 0 ) );
			spx = vertexList.get( 0 ).getPoints().get( 0 ).x;
			spy = vertexList.get( 0 ).getPoints().get( 0 ).y;
			spz = vertexList.get( 0 ).getPoints().get( 0 ).z;
			return 0;
		}

		//create empty adjacency and predecessor matrix

		/** the matrix that contains the length of the shortest path from vertex a to vertex b */
		double[][] adjacencyMatrix = new double[vertexList.size()][vertexList.size()];
		/** the matrix that contains the predecessor vertex of vertex b in the shortest path from vertex a to b */
		int[][] predecessorMatrix = new int[vertexList.size()][vertexList.size()];

		// initial conditions for both matrices
		/*
		 * create 2D-adjacency array with the distance between the nodes.
		 * distance i --> i = 0
		 * distance i -/> j = infinite (edge does not exist)
		 * distance i --> j = length of branch
		 */

		/*
		 * create 2D-predecessor array using interconnected vertices.
		 * the predecessor matrix contains the predecessor of j in a path from node i to j.
		 * distance i --> i = NIL (there is no edge between a single vertex)
		 * distance i -/> j = NIL (there is no edge between the vertices)
		 * distance i --> j = i (the initial matrix only contains paths with a single edge)
		 * 
		 * I'm using -1 as NIL since it cannot refer to an index (and thus a vertex)
		 */

		// applying initial conditions
		for(int i = 0 ; i < vertexList.size(); i++)
		{
			for(int j = 0 ; j < vertexList.size(); j++)
			{
				adjacencyMatrix[i][j]= Double.POSITIVE_INFINITY;
				predecessorMatrix[i][j]= -1;
			}
		}


		for (Edge edge : edgeList)
		{
			v1 = edge.getV1();
			v2 = edge.getV2();
			// use the index of the vertices as the index in the matrix
			row = vertexList.indexOf(v1);
			if(row == -1)
			{
				IJ.log("Vertex " + v1.getPoints().get(0) + " not found in the list of vertices!");
				continue;
			}

			column = vertexList.indexOf(v2);
			if(column == -1)
			{
				IJ.log("Vertex " + v2.getPoints().get(0) + " not found in the list of vertices!");
				continue;
			}

			/* 
			 * the diagonal is 0. 
			 * 
			 * Because not every vertex is a 'v1 vertex' 
			 * the [column][column] statement is needed as well.
			 *
			 * in an undirected graph the adjacencyMatrix is symmetric 
			 * thus A = Transpose(A)
			 */
			adjacencyMatrix[row][row] = 0;
			adjacencyMatrix[column][column] = 0;
			adjacencyMatrix[row][column] = edge.getLength();
			adjacencyMatrix[column][row] = edge.getLength();

			/* 
			 * the diagonal remains -1.
			 * for the rest I use the index of the vertex so I can later refer to the vertexList
			 * for the correct information
			 * 
			 * Determining what belongs where requires careful consideration of the definition
			 * 
			 * the array contains the predecessor of "column" in a path from "row" to "column"
			 * therefore in the other statement it is the other way around.   
			 */
			predecessorMatrix[row][row] = -1;
			predecessorMatrix[column][column] = -1;
			predecessorMatrix[row][column] = row;
			predecessorMatrix[column][row] = column;
		}
		// matrices now have their initial conditions


		// the warshall algorithm with k as candidate vertex and i and j walk through the adjacencyMatrix
		// the predecessor matrix is updated at the same time. 

		for(int k = 0 ; k < vertexList.size(); k++)
		{
			for(int i = 0 ; i < vertexList.size(); i++)
			{
				for(int j = 0 ; j < vertexList.size(); j++)
				{
					if(adjacencyMatrix[i][k] + adjacencyMatrix[k][j] < adjacencyMatrix[i][j])
					{
						adjacencyMatrix[i][j] = adjacencyMatrix[i][k] + adjacencyMatrix[k][j];
						predecessorMatrix[i][j] = predecessorMatrix[k][j];

					}
				}
			}
		}

		// find the maximum of all shortest paths
		for(int i = 0; i < vertexList.size(); i++)
		{
			for(int j = 0; j < vertexList.size(); j++)
			{
				// sometimes infinities still remain
				if (adjacencyMatrix[i][j] > maxPath && adjacencyMatrix[i][j] != Double.POSITIVE_INFINITY)
				{
					maxPath = adjacencyMatrix[i][j];
					a = i;
					b = j;

				}
			}
		}

		// trace back the longest shortest path
		reconstructPath(predecessorMatrix, a, b, edgeList, vertexList, shortestPathPoints);

		// !important return maxPath;
		return maxPath;

	}
	// end method warshallAlgorithm

	/**
	 * Reconstruction and visualisation of the longest shortest path found by the APSP warshall algorithm
	 *  
	 * @param predecessorMatrix the Matrix which contains the predecessor of vertex b in the shortest path from a to b
	 * @param startIndex the index of the row which contains the longest shortest path
	 * @param endIndex the index of the column which contains the longest shortest path
	 * @param edgeList the list of edges 
	 * @param vertexList the list of vertices
	 * @param shortestPathPoints contains points of the longest shortest path for each graph
	 * @author Huub Hovens
	 */
	private void reconstructPath(
			int[][] predecessorMatrix, 
			int startIndex, 
			int endIndex, 
			ArrayList<Edge> edgeList, 
			ArrayList<Vertex> vertexList,
			ArrayList<Point> shortestPathPoints)
	{
		// We know the first and last vertex of the longest shortest path, namely a and b
		// using the predecessor matrix we can now determine the path that is taken from a to b
		// remember a and b are indices and not the actual vertices.

		int b = endIndex;
		final int a = startIndex;

		while (b != a)
		{
			Vertex predecessor = vertexList.get(predecessorMatrix[a][b]);
			Vertex endvertex = vertexList.get(b);
			ArrayList< Edge > sp_edgeslist = new ArrayList< Edge >();
			Double lengthtest = Double.POSITIVE_INFINITY;
			Edge shortestedge = null;

			// search all edges for a combination of the two vertices
			for (Edge edge : edgeList)
			{

				if ((edge.getV1()==predecessor && edge.getV2()==endvertex) || (edge.getV1()==endvertex && edge.getV2()==predecessor))
				{
					// sometimes there are multiple edges between two vertices so add them to a list
					// for a second test
					sp_edgeslist.add(edge);
				}

			}
			// the second test
			// this test looks which edge has the shortest length in sp_edgeslist
			for (Edge edge : sp_edgeslist)
			{
				if (edge.getLength() < lengthtest)
				{
					shortestedge = edge;
					lengthtest = edge.getLength();
				}

			}
			// add vertex 1 points
			Vertex v1 = shortestedge.getV2() != predecessor ?
						shortestedge.getV2() : shortestedge.getV1();
			for (Point p : v1.getPoints())
			{
				if( ! shortestPathPoints.contains( p ))
				{
					shortestPathPoints.add(p);
					//setPixel(this.shortPathImage, p.x, p.y, p.z, SHORTEST_PATH);
				}
			}

			// add slab points of the shortest edge to the list of points
			ArrayList<Point> slabs = shortestedge.getSlabs();
			// reverse order if needed
			if( shortestedge.getV2() != predecessor )
				Collections.reverse( slabs );
			for (Point p : slabs)
			{
				shortestPathPoints.add(p);
				setPixel(this.shortPathImage, p.x, p.y, p.z, SHORTEST_PATH);
			}

			// add vertex 2 points too
			Vertex v2 = shortestedge.getV2() != predecessor ?
					shortestedge.getV1() : shortestedge.getV2();
			for (Point p : v2.getPoints())
			{
				if( ! shortestPathPoints.contains( p ))
				{
					shortestPathPoints.add(p);
					//setPixel(this.shortPathImage, p.x, p.y, p.z, SHORTEST_PATH);
				}
			}

			// now make the index of the endvertex the index of the predecessor so that the path now goes from
			// a to predecessor and repeat cycle
			b = predecessorMatrix[a][b];
		}
		if (shortestPathPoints.size() != 0)
		{
			this.spx = shortestPathPoints.get(0).x;
			this.spy = shortestPathPoints.get(0).y;
			this.spz = shortestPathPoints.get(0).z;
		}

	}
	// end method reconstructPath

	/**
	 * Get the stack containing all the trees labeld with their corresponding
	 * skeleton id.
	 *  
	 * @return labeled-skeleton image stack
	 */
	public ImageStack getLabeledSkeletons()
	{
		return labeledSkeletons;
	}

	private void createSettingsDialog() {
		settingsDialog = new GenericDialog("Analyze Skeleton");
		Font headerFont = new Font("SansSerif", Font.BOLD, 12);

		settingsDialog.setInsets(0, 0, 0);
		settingsDialog.addMessage("Elimination of Loops:", headerFont);
		settingsDialog.addChoice("Prune cycle method: ", AnalyzeSkeleton_.pruneCyclesModes,
				AnalyzeSkeleton_.pruneCyclesModes[pruneIndex]);

		settingsDialog.setInsets(20, 0, -15); //default top inset for 1st checkbox is 15
		settingsDialog.addMessage("Elimination of End-points:", headerFont);
		settingsDialog.addCheckbox("Prune ends", pruneEnds);
		settingsDialog.addCheckbox("Exclude ROI from pruning", protectRoi);

		settingsDialog.setInsets(20, 0, 0); //default top inset for subsequent checkboxes is 0
		settingsDialog.addMessage("Results and Output:", headerFont);
		settingsDialog.addCheckbox("Calculate largest shortest path", calculateShortestPath);
		settingsDialog.addCheckbox("Show detailed info", verbose);
		settingsDialog.addCheckbox("Display labeled skeletons", displaySkeletons);

		settingsDialog.addHelp(HELP_URL);
		dialogItemChanged(settingsDialog, null);
	}

	private void loadDialogSettings() {
		String index = Prefs.get(PRUNE_MODE_INDEX_KEY, String.valueOf(DEFAULT_PRUNE_MODE_INDEX));
		pruneIndex = Integer.parseInt(index); // fails to find the key
		pruneEnds = Prefs.get(PRUNE_ENDS_KEY, DEFAULT_PRUNE_ENDS);
		calculateShortestPath = Prefs.get(CALCULATE_PATH_KEY, DEFAULT_CALCULATE_SHORTEST_PATH);
		verbose = Prefs.get(VERBOSE_KEY, DEFAULT_VERBOSE);
		displaySkeletons = Prefs.get(DISPLAY_SKELETONS_KEY, DEFAULT_DISPLAY_SKELETONS);
	}

	private void saveDialogSettings() {
		Prefs.set(PRUNE_MODE_INDEX_KEY, String.valueOf(pruneIndex));
		Prefs.set(PRUNE_ENDS_KEY, pruneEnds);
		Prefs.set(CALCULATE_PATH_KEY, calculateShortestPath);
		Prefs.set(VERBOSE_KEY, verbose);
		Prefs.set(DISPLAY_SKELETONS_KEY, displaySkeletons);
	}

	private void setSettingsFromDialog() {
		pruneIndex = settingsDialog.getNextChoiceIndex();
		pruneEnds = settingsDialog.getNextBoolean();
		protectRoi = settingsDialog.getNextBoolean();
		calculateShortestPath = settingsDialog.getNextBoolean();
		verbose = settingsDialog.getNextBoolean();
		displaySkeletons = settingsDialog.getNextBoolean();
	}
}// end class AnalyzeSkeleton_
