package vnt;
/**
 * Node_Analysis.java - v1.0
 * Began: July 1, 2005
 * Last Updated: July 28, 2005 
 * 
 * Copyright (C) 2005 - Michael D. Miller - mdm162@truman.edu
 * Truman State University
 * 100 E. Normal
 * Kirksville, MO - 63501
 * 
 * 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 2
 * 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, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
 */

import ij.*;
import ij.gui.*;

import java.awt.*;
import java.awt.Color.*;
import ij.plugin.filter.PlugInFilter;
import ij.process.*;

import java.io.*;
import ij.io.*;
import java.lang.String.*;
import java.lang.Math.*;
import ij.plugin.filter.*;
import ij.measure.*; // needed for ResultsTable

import Coordinate.*;
import java.util.Stack;
import Node.*;
import Node_LinkedList.*;

//import java.io.File.*;
//import java.io.FileInputStream.*;
//import java.io.DataInputStream.*;

/**
 * <p>Multi-purpose class for final analysis.
 * Accepts a pruned, thinned distance map skeleton.
 * Performs the following operations:</p>
 * 
 * <ol>
 * <li>Generates node information at endpoints and branch points.</li>
 * <li>Traverses the thinned skeleton and computes adjacency list information.</li> 
 * <li>Merges nodes with significant overlap.</li>
 * <li>Outputs the final coordinate and adjacency list information to text files for Mathematica analysis.</li>
 * <li>Performs quantitative and topological analysis on the image's skeleton.</li>
 * </ol>
 * 
 * @author Michael Miller - Truman State University
 * @version 1.0
 * @since 1.0
 */
public class Node_Analysis extends VascularNetworkToolkit implements PlugInFilter {

    /**
     * An array of nodes. Holds all the nodes of the image. 
     */
    private Node [] nodes;

    /**
     * Our own persistent results table. This is responsible for displaying
     * the "Results" window data.
     */
    private static ResultsTable rt;
    
    /**
     * 
     */
    public Node_Analysis( ) {
        if (rt == null || !IJ.isResultsWindow())
            rt = new ResultsTable();
    }
    
    /**
    * Specifies the preconditions for the plug-in.
    * If this method succeeds then run() is called.
    *
    * <p>Pre: ImageJ is running and an 8-bit grayscale image is open. The plug-in was just activated.
    * <br />Post: Either an argument was processed, the image was not saved to a local folder, or the plug-in is cleared to run on the image.
    * @param arg Required by the interface. The argument list passed to the plug-in.
    * @param imp Required by the interface. The ImagePlus datastructure holding (among other things) information to grab path and filename.
    * @return If DONE is returned, ImageJ quits without run()'ing the plug-in. Otherwise, the plug-in signals to ImageJ that this plugin only handles 8-bit (256 grayscale) and will not change the original image.
    * @see #run
    */
    public int setup(String arg, ImagePlus ip) {
        if (arg.equals("about")) { 
            showAbout("Node Analysis","  * Calculates quantitative information from a given grayscale skeleton image. (Copyright 2005. Michael Miller mdm162@truman.edu)");
            return DONE; // exit without run()
        }

        // will only save if this succeeds
        getFileInformation(ip);
        
        return DOES_8G+NO_CHANGES; // success, run()
    }

    /**
     * Receives a binary image (specifically a thinned, pruned distance map
     * skeletonization) and computes information regarding the topological
     * properties of the image.
     * 1) Parses the image for possible node locations.
     * 2) Accesses the distance map to calculate node radii information.
     * 3) Simplifies the structures by testing for node overlap.
     * 4) Performs graph theoretic analysis based on the information generated. 
     *  
     * <p>Pre: The image was cleared to run by the setup() method.
     * <br />Post: Quantitative data is finally extracted from the image. 
     * @param bp Required by the interface. The access information to the original image.
     * @see #setup
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     * @see #mergeOverlappingNodes
     * @see #drawNodes
     * @see #drawSkeletonWithNodes
     * @see #drawNodeGraph
     * @see #writeTextFiles
     */
    public void run(ImageProcessor bp){ 
        // get the width and height information
        getDimensions(bp);
        // find all the nodes
        if (displayDebugText)
            System.out.print("About to generate all node information...\n");
        generateNodeInformation(bp);
        // find all the adjacencies
        if (displayDebugText)
            System.out.print("About to generate all adjacency information...\n");
        generateAdjacencyLists(bp);
        // display the nodes before merging
        //drawSkeletonWithNodes(bp, true);
        //drawSkeletonWithNodes(bp, false);
        // node counter before

        // display the nodes after merging
        if (displayDebugText)
            System.out.print("Drawing a node overlay...\n");
        drawNodes(false);
        // save the node display image
        saveFile("nodes");      

        // display the nodes graph
        //drawNodeGraph(bp,true);
        if (displayDebugText)
            System.out.print("Drawing a graph overlay...\n");
        drawNodeGraph(bp,false);
        // save the node/adjacencylist display image
        saveFile("graph");

        if (displayDebugText)
            System.out.print("About to merge overlapping nodes.\nBefore the merging: "+countNodes()+"\n");
        // collapse overlapping nodes onto themselves while preserving adjacency list information
        mergeOverlappingNodes(0);
        // node counter after
        if (displayDebugText)
            System.out.print("After the merging: "+countNodes()+"\n");
        if (displayDebugText)
            System.out.print("Drawing a new merged graph  and saving it..");
        drawNodeGraph(bp,false);
        saveFile("merged");

        // analyze results
        if (displayDebugText)
            System.out.print("Writing data to text files...\n");
        writeTextFiles();
        
        if (displayDebugText)
            System.out.print("Generating data to display to Results table...\n");
        displayProperties(bp);
        if (displayDebugText)
            System.out.print("Dumping text to the table...\n");
        displayResults();
        if (displayDebugText)
            System.out.print("...all done with Node_Analysis\n");
        return; 
    }
    
    /**
     * Parses the given image for nodes. Nodes are identifiable as pixels with
     * either one or three 8-neighbor backgorund 4-connected components.
     * Ie, nodes are either at endpoints or at branch points.  
     *  
     * <p>Pre: The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be WHITE-radius.
     * <br />Post: The local member nodes[] is created and populated with the nodes in the image.
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @return Returns the number of nodes that were found in the image.
     * @see #nodes
     */
    public int generateNodeInformation(ImageProcessor bp) {
        int x=0, y=0, count=0, components=0;
        Coordinate skeleton;
        // load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // get a count of the number of potential nodes
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                // only concerned about skeletal pixels
                if (pixel[x][y] != WHITE) {
                    skeleton = new Coordinate(x,y,pixel[x][y]);
                    components = getWhiteComponents(skeleton); 
                    // are there 1 (endpoint) or 3 (branch pt.) components?
                    if ((components == 1) || (components >= 3)) {
                        // node!
                        count ++;
                    }
                }
            }
        }
        // create a datastructure that holds them all
        nodes = new Node [count];
        // build the array of potential nodes
        count = 0;
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                // only concerned about skeletal pixels
                if (pixel[x][y] != WHITE) {
                    skeleton = new Coordinate(x,y,pixel[x][y]);
                    components = getWhiteComponents(skeleton); 
                    // are there 1 (endpoint) or 3 (branch pt.) components?
                    if ((components == 1) || (components >= 3)) {
                        // node!
                        // radius is the distance from white
                        nodes[count] = new Node(x,y,WHITE-pixel[x][y]);
                        count++; 
                    }
                }
            }
        }
        return count;
    }
    
    /**
     * Builds the adjacency list information by parsing the given image
     * along the skeleton. Each node's 8-neighbors are traversed until
     * another node is found.  
     *  
     * <p>Pre: The local member nodes[] is created and populated with the nodes in the image. The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be (WHITE-radius).
     * <br />Post: The nodes[] structure now contains all adjacency list information. 
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @see #generateNodeInformation
     */
    public void generateAdjacencyLists(ImageProcessor bp) {
        int x=0, y=0, i=0, j=0, color=0, components=0, total=0;
        Coordinate skeleton;
        // traverse the nodes and build adjacency list information
        // traverse the image and build node information
        Stack searchStack = new Stack();
        for(i=0; i<nodes.length; i++) {
            searchStack.clear();
//          searchStack.push(new Coordinate(nodes[i].getX(),nodes[i].getY(), nodes[i].getRadius()));
            setColor(nodes[i].getX(), nodes[i].getY(), BLACK);
            if (skeletalNeighbors(nodes[i].getX(),nodes[i].getY()) != 2) {
                for (x=nodes[i].getX()-1; x<=nodes[i].getX()+1; x++) {
                    for (y=nodes[i].getY()-1; y<=nodes[i].getY()+1; y++) {
                        color = getColor(x,y);          
                        if ((color != WHITE) && (color != BLACK)) {
                            searchStack.push(new Coordinate(x, y, color));
                        }
                    }
                }               
            } else { // a node is nearby if there are 2 at the start
                for (x=nodes[i].getX()-1; x<=nodes[i].getX()+1; x++) {
                    for (y=nodes[i].getY()-1; y<=nodes[i].getY()+1; y++) {
                        color = getColor(x,y);          
                        if ((color != WHITE) && (color != BLACK)) {
                            components = getWhiteComponents(new Coordinate(x,y,color));
                            if (((components == 1) || (components >= 3)) && (!( (nodes[i].getX()==x)&&(nodes[i].getY()==y) ) )) {
                                // node! adjacent!
                                for (j=0; j<nodes.length; j++) {
                                    if ((nodes[j].getX() == x) && (nodes[j].getY() == y)) {
                                        if (i != j) {
                                            nodes[i].addAdjacency(nodes[j]);
                                        }
                                        break; // j = count?
                                    }
                                }
                            }
                        }
                    }
                }
            }
            while (!searchStack.isEmpty()) {
                skeleton = (Coordinate)searchStack.pop();
                components = getWhiteComponents(skeleton); 
                x = skeleton.getX();
                y = skeleton.getY();
                // make sure we didn't already place this pixel
                if (pixel[x][y] != BLACK) {
                    // is this a node?
                    if (((components == 1) || (components >= 3)) && (!( (nodes[i].getX()==x)&&(nodes[i].getY()==y) ) )) {
                        // node! adjacent!
                        // search all nodes and find the one we hit
                        for (j=0; j<nodes.length; j++) {
                            if ((nodes[j].getX() == x) && (nodes[j].getY() == y)) {
                                nodes[i].addAdjacency(nodes[j]);
                                break; // j = count?
                            }
                        }
                    } else if (skeletalNeighbors(x,y) != 2) {
                        // about to hit a node, find it and quit
                        pixel[x][y] = BLACK;
                        total = 0;
                        for (x=skeleton.getX()-1; x<=skeleton.getX()+1; x++) {
                            for (y=skeleton.getY()-1; y<=skeleton.getY()+1; y++) {
                                color = getColor(x,y);          
                                if ((color != WHITE) && (color != BLACK)) {
                                    components = getWhiteComponents(new Coordinate(x,y,color));
                                    if (((components == 1) || (components >= 3)) && (!( (nodes[i].getX()==x)&&(nodes[i].getY()==y) ) )) {
                                        // node! adjacent!
                                        for (j=0; j<nodes.length; j++) {
                                            if ((nodes[j].getX() == x) && (nodes[j].getY() == y)) {
                                                if (i != j) {
                                                    nodes[i].addAdjacency(nodes[j]);
                                                    total ++;
                                                }
                                                break; // j = count?
                                            }
                                        }
                                    }
                                }
                            }
                        }
                        if (total == 0) { // found none, must have detected we were close to a starter 3-way                            
                            pixel[skeleton.getX()][skeleton.getY()] = BLACK;
                            for (x=skeleton.getX()-1; x<=skeleton.getX()+1; x++) {
                                for (y=skeleton.getY()-1; y<=skeleton.getY()+1; y++) {
                                    color = getColor(x,y);          
                                    if ((color != WHITE) && (color != BLACK)) {
                                        searchStack.push(new Coordinate(x, y, color));
                                    }
                                }
                            }
                        }
                    } else {
                        // not a node, keep looking
                        pixel[x][y] = BLACK;
                        // check the cardinal directions to add 4-connected
                        // neighbors that satisfy the high threshold
                        
                        for (x=skeleton.getX()-1; x<=skeleton.getX()+1; x++) {
                            for (y=skeleton.getY()-1; y<=skeleton.getY()+1; y++) {
                                color = getColor(x,y);          
                                if ((color != WHITE) && (color != BLACK)) {
                                    searchStack.push(new Coordinate(x, y, color));
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    /**
     * Counts the number of nodes in the image.  
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     */
    public int countNodes() {
        int total=0, i;
        
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                total ++;
            }
        }
        return total;
    }

    /**
     * Counts the number of edges in the image.   
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated. It helps a lot if the generateAdjacencyLists procedure has also been called.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public int countEdges() {
        int total=0, i;
        int x, y, adjacentX, adjacentY;
        Node_LinkedList adjacencies;
        
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                x = nodes[i].getX();
                y = nodes[i].getY();
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }

                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        adjacentX = adjacencies.getNode().getX();
                        adjacentY = adjacencies.getNode().getY();
                        // don't count reflexive adjacencies
                        if (!((x == adjacentX) && (y == adjacentY)))
                            total ++;
                        adjacencies = adjacencies.getAfter();
                    }
                }           
            }
        }
        total /= 2; // all edges are counted twice
        return total;
    }

    /**
     * Counts the number of branching points from the nodes[]
     * datastructure. A branching point is a node with two
     * or more adjacencies.   
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated. It helps a lot if the generateAdjacencyLists procedure has also been called.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public int countBranchingPoints() {
        int total=0, i=0, degrees=0, count=0;
        int x, y, adjacentX, adjacentY;
        Coordinate skeleton;
        Node_LinkedList adjacencies;

        // get a count of the number of potential nodes
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                x = nodes[i].getX();
                y = nodes[i].getY();
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }
                degrees = 0;
                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        adjacentX = adjacencies.getNode().getX();
                        adjacentY = adjacencies.getNode().getY();
                        // don't count reflexive adjacencies
                        if (!((x == adjacentX) && (y == adjacentY)))
                            degrees ++;
                        adjacencies = adjacencies.getAfter();
                    }
                }
                if (degrees >= 2) {
                    // branching point!
                    count ++;
                }
            }
        }
        return count;
    }
    
    /**
     * Counts the number of end points from the nodes[]
     * datastructure. An end point is a node with a single
     * adjacency.
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated. It helps a lot if the generateAdjacencyLists procedure has also been called.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public int countEndPoints() {
        int total=0, i=0, degrees=0, count=0;
        int x, y, adjacentX, adjacentY;
        Coordinate skeleton;
        Node_LinkedList adjacencies;

        // get a count of the number of potential nodes
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                x = nodes[i].getX();
                y = nodes[i].getY();
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }
                degrees = 0;
                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        adjacentX = adjacencies.getNode().getX();
                        adjacentY = adjacencies.getNode().getY();
                        // don't count reflexive adjacencies
                        if (!((x == adjacentX) && (y == adjacentY)))
                            degrees ++;
                        adjacencies = adjacencies.getAfter();
                    }
                }
                if (degrees == 1) {
                    // end point!
                    count ++;
                }
            }
        }
        return count;
    }
    
    /**
     * Counts the number of isolated points from the nodes[]
     * datastructure. An isolated point is a node with no
     * adjacencies.
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated. It helps a lot if the generateAdjacencyLists procedure has also been called.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public int countIsolatedPoints() {
        int total=0, i=0, degrees=0, count=0;
        int x, y, adjacentX, adjacentY;
        Coordinate skeleton;
        Node_LinkedList adjacencies;

        // get a count of the number of potential nodes
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                x = nodes[i].getX();
                y = nodes[i].getY();
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }
                degrees = 0;
                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        adjacentX = adjacencies.getNode().getX();
                        adjacentY = adjacencies.getNode().getY();
                        // don't count reflexive adjacencies
                        if (!((x == adjacentX) && (y == adjacentY)))
                            degrees ++;
                        adjacencies = adjacencies.getAfter();
                    }
                }
                if (degrees == 0) {
                    // isolated point!
                    count ++;
                }
            }
        }
        return count;
    }

    /**
     * Counts the number of edges in the image.
     * Returns the maximum degree.   
     *  
     * <p>Pre: The local member nodes[] is checked for the number of non-null nodes. The nodes[] array must be instantiated. It helps a lot if the generateAdjacencyLists procedure has also been called.
     * <br />Post: The nodes[] structure is searched for non-null nodes and returns the total count. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public int countMaxDegree() {
        int total=0, i, max=0;
        int x, y, adjacentX, adjacentY;
        Node_LinkedList adjacencies;
        
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                x = nodes[i].getX();
                y = nodes[i].getY();
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }
                total = 0;
                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        adjacentX = adjacencies.getNode().getX();
                        adjacentY = adjacencies.getNode().getY();
                        // don't count reflexive adjacencies
                        if (!((x == adjacentX) && (y == adjacentY)))
                            total ++;
                        adjacencies = adjacencies.getAfter();
                    }
                }   
                if (total > max) {
                    max = total;        
                }
            }
        }
        return max;
    }

    /**
     * Parses the given image for skeletal structure.  
     * Returns the total area in pixels.
     *  
     * <p>Pre: The ImageProcessor bp must exist. This should be the segmented image.
     * <br />Post: The local member nodes[] is created and populated with the nodes in the image.
     * @param bp The ImageProcessor for the segmented image.
     * @return Returns the total area of the vascular network in pixels.
     */
    public double calculateTotalNetworkArea(ImageProcessor bp) {
        int x=0, y=0;
        double area=0;
        double realWorldUnitsPixelArea = imageAspectRatio*distanceInUnitsPerPixel*distanceInUnitsPerPixel;
        // load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // get a count of the number of potential nodes
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                // only concerned about skeletal pixels
                if (getColor(x,y) != WHITE)
                    area += realWorldUnitsPixelArea;
            }
        }
        return area;
    }

    /**
     * Parses the given image for skeletal structure.  
     * Length of the network is given by the equation:
     * L = N + [sqrt(2)-1]*Nd - 1
     * Where N is the number of pixels of the skeleton
     * and Nd is the number of diagonally connected pairs of pixels.
     *  
     * <p>Pre: The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be WHITE-radius.
     * <br />Post: The local member nodes[] is created and populated with the nodes in the image.
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @return Returns the total network length of the skeleton.
     */
    public double calculateTotalNetworkLength(ImageProcessor bp) {
        int x=0, y=0;
        double L=0, N=0, Nd=0, Nv=0, Nh=0, Nmod;
        double pixelWidth,pixelHeight,pixelDiagonal;
        pixelWidth = 1;
        pixelHeight = imageAspectRatio;
        pixelDiagonal = Math.sqrt(pixelWidth+pixelHeight);
        // load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // get a count of the number of potential nodes
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                // only concerned about skeletal pixels
                if (getColor(x,y) != WHITE) {
                    N ++;
                    if (getColor(x-1,y) != WHITE)
                        Nh ++;
                    if (getColor(x+1,y) != WHITE)
                        Nh ++;
                    if (getColor(x,y-1) != WHITE)
                        Nv ++;
                    if (getColor(x,y+1) != WHITE)
                        Nv ++;
                    // now count the diagonally connected pairs
                    if (getColor(x-1,y-1) != WHITE) { // upper left
                        if (getColor(x-1,y) == WHITE && getColor(x,y-1) == WHITE)
                            Nd ++;
                    }
                    if (getColor(x+1,y+1) != WHITE) { // lower right
                        if (getColor(x+1,y) == WHITE && getColor(x,y+1) == WHITE)
                            Nd ++;
                    }
                    if (getColor(x-1,y+1) != WHITE) { // lower left
                    if (getColor(x-1,y) == WHITE && getColor(x,y+1) == WHITE)
                            Nd ++;
                    }
                    if (getColor(x+1,y-1) != WHITE) { // upper right
                    if (getColor(x+1,y) == WHITE && getColor(x,y-1) == WHITE)
                            Nd ++;
                    }
                }
            }
        }
        Nd /= 2; // counted each pair twice
        Nv /= 2; // counted each pair twice
        Nh /= 2; // counted each pair twice
        Nmod = N - Nd - Nv - Nh;
        // Note: Nmod should be the number of pixels.
        
        //L = N*imageAspectRatio + (Math.sqrt(1+imageAspectRatio)-1)*Nd - 1;
        //pixelWidth,pixelHeight,pixelDiagonal;
        L = Nd*pixelDiagonal + Nh*pixelWidth + Nv*pixelHeight;
        L *= distanceInUnitsPerPixel;
        return L;
    }

    /**
     * Takes the segmented image and counts the number of enclosed meshes.
     * This is largely taken from the FindEdges_Segment.java algorithm,
     * although this is optimized for counting.
     *  
     * <p>Pre: The image must have been converted to binary.
     * <br />Post: As long as a valid ImageProcessor is given, the number of enclosed meshes is returned. 
     * @param bp The ImageProcessor for the segmented image.
     * @return Returns -1 on error. Otherwise returns a non-negative number that represents the number of 4-connected enclosed meshes.
     */
    public int calculateMeshCount(ImageProcessor bp) {
        if (bp == null) {
            return -1;
        }
        // at this point assume that any meshes we find are of appropriate size
        int x=0, y=0, seekX=0, seekY=0, meshCount=0;
        boolean enclosed = true;
        Stack searchStack = new Stack();
        Coordinate seeker = new Coordinate(x,y,WHITE);
        // load our original image into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // begin searching, any meshes we find count
        for (y=0; y<height; y++) {
            for (x=0; x<width; x++) {
                if (getColor(x,y) == WHITE) {
                    searchStack.clear();
                    seeker = new Coordinate(x,y,WHITE);
                    searchStack.push(seeker);
                    enclosed = true;
                    while (!searchStack.isEmpty()) {
                        seeker = (Coordinate)searchStack.pop();
                        seekX = seeker.getX();
                        seekY = seeker.getY();
                        // make sure we didn't already place this pixel
                        if (getColor(seekX,seekY) == WHITE) {
                            setColor(seekX, seekY, GRAY);
                            // set up 4 connected search left and right
                            seekY = seeker.getY();
                            for (seekX=seeker.getX()-1; seekX<=seeker.getX()+1; seekX+=2) {
                                if (seekX<0 || seekX>=width) {
                                    enclosed = false;
                                } else {
                                    if (getColor(seekX,seekY) == WHITE) {
                                        searchStack.push(new Coordinate(seekX, seekY, WHITE));
                                    }
                                }
                            }
                            // set up 4 connected search up and down
                            seekX = seeker.getX();
                            for (seekY=seeker.getY()-1; seekY<=seeker.getY()+1; seekY+=2) {
                                if (seekY<0 || seekY>=height) {
                                    enclosed = false;
                                } else {
                                    if (getColor(seekX,seekY) == WHITE) {
                                        searchStack.push(new Coordinate(seekX, seekY, WHITE));
                                    }
                                }
                            }
                        }
                    }
                    if (enclosed) {
                        meshCount ++;
                    }
                }
            }
        }
        
        return meshCount;
    }   
    
    /**
     * Parses the given image for skeletal structure.
     * Stays within the given 8-connected component.  
     * Length of the network is given by the equation:
     * L = N + [sqrt(2)-1]*Nd - 1
     * Where N is the number of pixels of the skeleton
     * and Nd is the number of diagonally connected pairs of pixels.
     *  
     * <p>Pre: The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be WHITE-radius.
     * <br />Post: The local member nodes[] is created and populated with the nodes in the image.
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @return Returns the total network length of the skeleton.
     */
    public double calculateConnectedComponentNetworkLength(ImageProcessor bp, Coordinate startingPoint) {
        int x=0, y=0, color;
        Coordinate search;
        double L=0, N=0, Nd=0;
        // load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // get a count of the number of potential nodes

        Stack searchStack = new Stack();
        searchStack.clear();
        if (startingPoint.getColor() != WHITE)
            searchStack.push(startingPoint);
        while (!searchStack.isEmpty()) {
            search = (Coordinate)searchStack.pop();
            x = search.getX();
            y = search.getY();
            // make sure we didn't already place this pixel
            if (getColor(x, y) != BLACK) { // new pixel
                setColor(x, y, BLACK);
                N ++;
                // now count the diagonally connected pairs
                if (getColor(x-1,y-1) != WHITE) // upper left
                    Nd ++;
                if (getColor(x+1,y+1) != WHITE) // lower right
                    Nd ++;
                if (getColor(x-1,y+1) != WHITE) // lower left
                    Nd ++;
                if (getColor(x+1,y-1) != WHITE) // upper right
                    Nd ++;

                for (x=search.getX()-1; x<=search.getX()+1; x++) {
                    for (y=search.getY()-1; y<=search.getY()+1; y++) {
                        color = getColor(x,y);          
                        if ((color != WHITE) && (color != BLACK)) {
                            searchStack.push(new Coordinate(x, y, color));
                        }
                    }
                }               
            }
        }
        Nd /= 2; // counted each pair twice
        L = N + (Math.sqrt(2)-1)*Nd - 1;
        return L;
    }

    /**
     * Parses the given image for skeletal structure.
     * Returns the average skeletal width.
     * Note: This method does NOT account for aspect ratio. 
     *  
     * <p>Pre: The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be WHITE-radius.
     * <br />Post: The local member nodes[] is created and populated with the nodes in the image.
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @return Returns the average width of the skeleton.
     * @see #WHITE
     */
    public double calculateAverageSkeletalWidth(ImageProcessor bp) {
        int x=0, y=0;
        double average=0, total=0, thickness=0, count=0;
        // load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);
        // get a count of width of the skeleton
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                // only concerned about skeletal pixels
                thickness = WHITE - getColor(x,y); // assumes WHITE=255, the max greyscale color
                if (thickness > 0) {
                    total += thickness;
                    count ++;
                }
            }
        }
        average = distanceInUnitsPerPixel*total/count;
        return average;
    }

    /**
     * Performs node merging accross the entire image.
     * 
     * <p>Pre: The local member nodes[] is created and populated with the nodes in the image. Some nodes may be null.
     * <br />Post: The nodes[] structure now contains all adjacency list information. 
     * @param overlap The amount of overlap required before two nodes are considered the same. The parameter should be in the range of [0, 1], where 0 implies two tangent circles will merged, and 1 requires that a node be completely inside another node to merge. 
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public void mergeOverlappingNodes(double overlap) {
        int x=0, y=0, i=0, j=0, size=0;
        
        for (i=0; i<nodes.length; i++) {
            for (j=0; j<nodes.length; j++) {
                if ((i != j) && (nodes[i] != null) && (nodes[j] != null)) {
                    // squared x distance
                    x = (nodes[i].getX()-nodes[j].getX())*(nodes[i].getX()-nodes[j].getX());
                    // squared y distance
                    y = (nodes[i].getY()-nodes[j].getY())*(nodes[i].getY()-nodes[j].getY());
                    // squared circle radii distance
                    size = (nodes[i].getRadius()+nodes[j].getRadius())*(nodes[i].getRadius()+nodes[j].getRadius());
                    if ((x+y) < (1-overlap)*size) {
                        // overlap!
                        nodes[i].mergeNodes(nodes[j]);
                        //color = nodes[j].getAdjacencySize();
                        //components = nodes[j].dissolve();
                        nodes[j] = null;
                    }
                }
            }
        }

        // dissolve tiny nodes
/*
        for (i=0; i<nodes.length; i++) {
            if (nodes[i] != null) {
                if (nodes[i].getRadius() < 4) {
                    nodes[i].dissolve();
                    nodes[i] = null;
                }       
            }
        }
        */
    }
    
    /**
     * Displays the computed nodes.  
     *  
     * <p>Pre: The local member nodes[] is created and populated with the nodes in the image. Some nodes may be null (due to merging).
     * <br />Post: A new image should be created with circular nodes overlayed on a white background. 
     * @param drawSimple If true, the image will be drawn with constant sized nodes. If false, the image will be drawn with true node sizes.
     * @return Returns true on successful draw. False otherwise.
     * @see #generateNodeInformation
     */
    public boolean drawNodes(boolean drawSimple) {
        String title;
        if (drawSimple) {
            title = "Simple Nodes";
        }else {
            title = "Nodes";
        }
        ImagePlus imp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor nip = imp.getProcessor();
        // draws circles around nodes
        for (int i=0; i<nodes.length; i++) {
            // specify a circle at each node            
            if (nodes[i] != null) {
                if (drawSimple) {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),5,BLACK); 
                } else {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),nodes[i].getRadius(),BLACK); 
                }
            }
        }
        imp.show();
        IJ.selectWindow(title);
        return true;
    }

    /**
     * Displays the image with nodes overlayed.  
     *  
     * <p>Pre: The local member nodes[] is created and populated with the nodes in the image. Some nodes may be null (due to merging). The ImageProcessor bp must exist. This should be the pruned, thinned distance map skeleton. Background should be white, and the color of the skeleton should be (WHITE-radius).
     * <br />Post: A new image should be created of the original with circular nodes overlayed. 
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @param drawSimple If true, the image will be drawn with a black skeleton and constant sized nodes. If false, the image will be drawn with distance map skeleton and true node sizes.
     * @return Returns true on successful draw. False otherwise.
     * @see #generateNodeInformation
     */
    public boolean drawSkeletonWithNodes(ImageProcessor bp, boolean drawSimple) {
        int x=0, y=0, i=0;
        String title;
        // re-load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);

        if (drawSimple) {
            title = "Simple Skeleton and Nodes";
        } else {
            title = "Skeleton with Nodes";
        }
        ImagePlus imp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor nip = imp.getProcessor();
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                if (drawSimple) {
                    if (pixel[x][y] != WHITE) {
                        nip.putPixel(x,y,BLACK);
                    }
                } else {
                    nip.putPixel(x,y,pixel[x][y]);
                }
            }
        }
        // draws circles around nodes
        int size;
        for (i=0; i<nodes.length; i++) {
            // specify a circle at each node            
            if (nodes[i] != null) {
                if (drawSimple) {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),5,GRAY); 
                } else {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),nodes[i].getRadius(),GRAY); 
                }
            }
        }
        imp.show();
        IJ.selectWindow(title);
        return true;
    }

    /**
     * Displays the graph of the with nodes with their edges.  
     *  
     * <p>Pre: The local member nodes[] is created and populated with the nodes in the image. Some nodes may be null (due to merging).
     * <br />Post: A new image should be created of circular nodes drawn with their edges.
     * @param bp The ImageProcessor for the pruned, thinned distance map skeleton.
     * @param drawSimple If true, the image will be drawn with constant sized nodes. If false, the image will be drawn with true node sizes.
     * @return Returns true on successful draw. False otherwise.
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public boolean drawNodeGraph(ImageProcessor bp, boolean drawSimple) {
        int x=0, y=0, targetX, targetY, i=0;
        String title;
        Node_LinkedList adjacencies;
        Line myLine;
        // re-load the distance-map skeleton into memory
        pixel = VascularNetworkToolkit.LoadImage(bp);

        if (drawSimple) {
            title = "Simple Graph";
        } else {
            title = "Graph";        
        }
        ImagePlus imp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor nip = imp.getProcessor();

        // draws circles around nodes
        int size;
        for (i=0; i<nodes.length; i++) {
            // specify a circle at each node            
            if (nodes[i] != null) {
                if (drawSimple) {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),5,GRAY); 
                } else {
                    VascularNetworkToolkit.drawCircle(nip,nodes[i].getX(),nodes[i].getY(),nodes[i].getRadius(),GRAY); 
                }
                adjacencies = nodes[i].getAdjacency();
                while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                    adjacencies = adjacencies.getBefore();
                }

                while (adjacencies != null) {
                    if (adjacencies.getNode() != null) {
                        targetX = adjacencies.getNode().getX();
                        targetY = adjacencies.getNode().getY();
                        VascularNetworkToolkit.drawLine(nip, nodes[i].getX(), nodes[i].getY(), targetX, targetY, BLACK);
                        adjacencies = adjacencies.getAfter();
                    }
                }           
            }
        }
        imp.show();
        IJ.selectWindow(title);
        return true;
    }

    /**
     * Simulates the original code by Jason Hart and Chris Miller.
     * Dumps coordinates and adjacency list information to a text file
     * for analysis in Mathematica using Kathleen Field's code.
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public void writeTextFiles () {
        FileOutputStream out; // declare a file output object
        PrintStream p; // declare a print stream object
        int i, j, f, nodeCount;
        int [] nodeListing;
        Node adjacent;
        Node_LinkedList adjacencies;
        boolean firstEntry, veryFirstEntry;
        // i is the memory index
        // here the first node is often 0, some may be null
        // f is the actual number associated +1
        // so the first node is 1, the second is 2, etc
        
        // create coord-prefix file
        try{

            String fileName = new String(originalName.toCharArray(),0,originalName.length()-4);
            String pathName = directory + "coord"+fileName+".txt";
            // Create a new file output stream
            out = new FileOutputStream(pathName);
            // Connect print stream to the output stream
            p = new PrintStream( out );
            // NOTE: build an array of i -> f conversions
            nodeCount = 0;
            for (i=0; i<nodes.length; i++) {
                if (nodes[i] != null) {
                    nodeCount ++;
                }
            }
            nodeListing = new int [nodes.length];
            nodeCount = 0;
            for (i=0; i<nodes.length; i++) {
                // specify a circle at each node
                nodeListing[i] = -1;            
                if (nodes[i] != null) {
                    nodeCount ++;
                    nodeListing[i] = nodeCount;
                }
            }
            p.print(nodeCount+"\n");// print the number of nodes
            f = 0;
            for(i=0;i<nodes.length;i++) {
                if (nodes[i] != null) {
                    //For each node print the xy coordinates and number of adjacencies.
                    p.print(nodes[i].getX()+" "+nodes[i].getY()+" "+(nodes[i].getAdjacencySize())+"\n");
                    // display the adjacency list
                    
                    // get the first
                    adjacencies = nodes[i].getAdjacency();
                    while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                        adjacencies = adjacencies.getBefore();
                    }
                    // iterate through and display them all
                    while (adjacencies != null) {
                        adjacent = adjacencies.getNode();
                        if (adjacent != null) {
                            // note: this is terrible!!! fix this!!
                            for (j=0; j<nodes.length; j++) {
                                if (adjacent == nodes[j]) {
                                    p.print(" " + nodeListing[j]);
                                }
                            }
                            // ::end note
                        }
                        adjacencies = adjacencies.getAfter();
                    }           
                    p.print("\n");
                }
            }
            p.close();
            out.close();
        }catch (Exception e) {
            System.out.print("Error writing to file\n");
        }

        // create skel-prefix file      
        try{

            String fileName = new String(originalName.toCharArray(),0,originalName.length()-4);
            String pathName = directory + "skel"+fileName+".txt";
            // Create a new file output stream
            out = new FileOutputStream(pathName);
            // Connect print stream to the output stream
            p = new PrintStream( out );
            // NOTE: build an array of i -> f conversions
            nodeCount = 0;
            for (i=0; i<nodes.length; i++) {
                if (nodes[i] != null) {
                    nodeCount ++;
                }
            }
            nodeListing = new int [nodes.length];
            nodeCount = 0;
            for (i=0; i<nodes.length; i++) {
                // specify a circle at each node
                nodeListing[i] = -1;            
                if (nodes[i] != null) {
                    nodeCount ++;
                    nodeListing[i] = nodeCount;
                }
            }
            p.print("{");
            f = 0;
            veryFirstEntry = true;
            for(i=0;i<nodes.length;i++) {
                if (nodes[i] != null) {
                    //For each node print the xy coordinates and number of adjacencies.
                    if (veryFirstEntry) {
                        p.print("{");
                        veryFirstEntry = false;
                    } else {
                        p.print(",{");
                    }
                    // display the adjacency list
                    
                    // get the first
                    firstEntry = true;
                    adjacencies = nodes[i].getAdjacency();
                    while ((adjacencies != null) && (adjacencies.getBefore() != null)) {
                        adjacencies = adjacencies.getBefore();
                    }
                    // iterate through and display them all
                    while (adjacencies != null) {
                        adjacent = adjacencies.getNode();
                        if (adjacent != null) {
                            // note: this is terrible!!! fix this!!
                            for (j=0; j<nodes.length; j++) {
                                if (i != j) {
                                    if (adjacent == nodes[j]) {
                                        if (firstEntry) {
                                            firstEntry = false;
                                            p.print(nodeListing[j]);
                                        } else {
                                            p.print("," + nodeListing[j]);
                                        }
                                    }
                                }
                            }
                            // ::end note
                        }
                        adjacencies = adjacencies.getAfter();
                    }           
                    p.print("}");
                }
            }
            p.print("}");
            p.close();
            out.close();
        }catch (Exception e) {
            System.out.print("Error writing to file\n");
        }       
    }
    
    /**
     * Dumps all quantitative and topological properties extractable
     * from the image into the ImageJ results table.
     * @see #generateNodeInformation
     * @see #generateAdjacencyLists
     */
    public ResultsTable displayProperties (ImageProcessor skeleton) {
        rt.incrementCounter();
        int column = 1;

        //quantitative properties
        IJ.open(directory+"segmented_"+originalName);
        ImagePlus segmentedImage = IJ.getImage();
        ImageProcessor segmentedIP = segmentedImage.getProcessor();
        double totalArea = calculateTotalNetworkArea(segmentedIP);
        int meshCount = calculateMeshCount(segmentedIP);
        IJ.run("Close");
        double totalLength = calculateTotalNetworkLength(skeleton);
        double averageWidth = calculateAverageSkeletalWidth(skeleton);
        double minMeshArea = maximumAreaToIgnoreMesh;
        minMeshArea *= distanceInUnitsPerPixel*distanceInUnitsPerPixel;
        rt.addLabel("File Name + Unit", originalName + " + " + distanceUnitName);
        rt.addValue(column, distanceInUnitsPerPixel);
        rt.setHeading(column++,"Units per Pixel");
        rt.addValue(column, imageAspectRatio);
        rt.setHeading(column++,"Pixel Aspect Ratio");
        rt.addValue(column, totalArea);
        rt.setHeading(column++,"Vascular Area");
        rt.addValue(column, averageWidth);
        rt.setHeading(column++,"Average Vessel Width");
        rt.addValue(column, totalLength);
        rt.setHeading(column++,"Network Length");
        rt.addValue(column, minMeshArea);
        rt.setHeading(column++,"Mesh Minimum Area");
        rt.addValue(column, meshCount);
        rt.setHeading(column++,"Enclosed Mesh Count");
        
        // topological properties
        int nodeCount = countNodes();
        int branchingPoints = countBranchingPoints();
        int endPoints = countEndPoints();
        int isolatedPoints = countIsolatedPoints();
        int edgeCount = countEdges();
        int highestDegree = countMaxDegree();
        
        rt.addValue(column, nodeCount);
        rt.setHeading(column++,"Total Node Count");
        rt.addValue(column, branchingPoints);
        rt.setHeading(column++,"Branching Points");
        rt.addValue(column, endPoints);
        rt.setHeading(column++,"End Points");
        rt.addValue(column, isolatedPoints);
        rt.setHeading(column++,"Isolated Points");
        rt.addValue(column, edgeCount);
        rt.setHeading(column++,"Edge Count");
        rt.addValue(column, highestDegree);
        rt.setHeading(column++,"Highest Degree");
        return rt;
    }

    /** 
     * Writes the last row in the results table to the ImageJ window.
     * (Parts of code stolen from ImageJ's Analyzer.java) 
     **/
    public void displayResults() {
        int counter = rt.getCounter();
        if (!IJ.isResultsWindow()) {
            IJ.setColumnHeadings(rt.getColumnHeadings());
        }
        String resultsLine = rt.getRowAsString(counter-1);  
        // display to the ImageJ table
        IJ.write(resultsLine);
        // write to a local file
        FileOutputStream out; // declare a file output object
        PrintStream p; // declare a print stream object
        String fileName;
        String pathName;
        try{
            fileName = new String(originalName.toCharArray(),0,originalName.length()-4);
            pathName = directory +"results_"+fileName+".txt";
            // Create a new file output stream
            out = new FileOutputStream(pathName);
            // Connect print stream to the output stream
            p = new PrintStream( out );
            // print the headings and results
            p.println(rt.getColumnHeadings());
            p.println(resultsLine); 
            p.close();
            out.close();
        }catch (Exception e) {
            System.out.print("Error writing to file\n");
        }       

        // test if the results file already exists
        fileName = "results.txt";
        pathName = directory + fileName;
        boolean resultsExist = true;        
        try {
            FileInputStream testInput = new FileInputStream(pathName);
            // didn't crash, file already exists
            testInput.close();
        } catch (Exception e) {
            // the file doesn't exist!
            resultsExist = false;
        }
        // now write to it
        try{
            // Create a new file output stream to append
            out = new FileOutputStream(pathName, resultsExist);
            // Connect print stream to the output stream
            p = new PrintStream( out );
            // if the file is empty, display the column headings
            if (!resultsExist) {
                p.println(rt.getColumnHeadings());
            }
            // print the results
            p.println(resultsLine); 
            p.close();
            out.close();
            
        }catch (Exception e) {
            System.out.print("Error writing to file\n");
        }       
    }
}