package vnt;
/**
 * DistanceMap_Skeleton.java - v1.0
 * Began: June 29, 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.*;

// Voronoi/Delaunay code required stuff
import java.util.Iterator;
import Delaunay.DelaunayTriangulation;
import Delaunay.Pnt;
import Delaunay.Simplex;
// my stuff
import Coordinate.*;
import java.util.Stack; // needed for searching

/**
 * <p>Uses ImageJ to generate a Distance Map for a segmented image.
 * Attempts to trace the darkest ridges of the distance map, providing
 * a closer approximation to the Blum Medial Axis Transform than the
 * Binary -> Skeletonization method (a thinning algorithm) given by
 * ImageJ.</p>
 * 
 * <p>The computed skeleton is not binary, it retains the distance map
 * information along the skeleton. This is (in theory) a one-to-one
 * invertible transform from the binary segmented image to the skeleton.</p>
 * 
 * @author Michael Miller - Truman State University
 * @version 1.0
 * @since 1.0
 */
public class DistanceMap_Skeleton extends VascularNetworkToolkit implements PlugInFilter {

    /**
    * 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("Distance Map Skeletonization","  * Computes the distance map from a given binary segmented image, and traces the ridges to compute a grayscale skeleton. (Copyright 2005. Michael Miller mdm162@truman.edu)");
            return DONE; // exit without run()
        }

        // will only save if this succeeds
        getFileInformation(ip);

        if (ip != null) {
            ImageStatistics stats=ip.getStatistics();
            if (stats != null) {
                if (stats.histogram[0]+stats.histogram[255]!=stats.pixelCount){
                    IJ.error("8-bit binary image (0 and 255) required.");
                    return DONE; // exit without run()
                }
            }
        }

        return DOES_8G+NO_CHANGES; // success, run()
    }

    /**
     * Receives a binary image and computes an approximation of the
     * Blum Medial Axis Transform on it.
     * 1) Performs distance map operation on the segementation image.
     * 2) Trace the ridges of the distance map.
     * 3) Generates a skeleton of the ridges of the distance map (which retain the color information from the distance map).
     *  
     * <p>Pre: The image was cleared to run by the setup() method.
     * <br />Post: The image is processed by the distance map and ridge tracing routines. A new distance map skeleton image is drawn. 
     * @param bp Required by the interface. The access information to the original image.
     * @see #setup
     * @see #generateDistanceMap
     * @see #generateRidgeTrace
     * @see #generateThinnedSkeleton
     */
    public void run(ImageProcessor bp){ 
        // get the width and height information
        getDimensions(bp);
        // generate the distance map
        ImageProcessor distanceMap = generateDistanceMap(bp);
        
        // Algorithm 1 - Produces a disconnected skeleton
        // Note: Insightful, but produces incomplete skeleton.
        // Attempts to trace along the ridges of the distance map.
        // this sadly doesn't work despite my best efforts ;-(
        //ImageProcessor ridgeTrace = generateRidgeTrace(distanceMap);
        // thin the grayscale skeleton using the Binary -> Skeletonization routine
        // ImageProcessor finalSkeleton = generateThinnedSkeleton(ridgeTrace);
        
        // Algorithm 2 - The standard "thinning" algorithm that ImageJ uses. 
        ImageProcessor cheatingSkeleton = generateThinningSkeletonWithEDM(distanceMap);
        
        // Algorithm 3 - Implementation of the Voronoi/Delaunay boundaries
        //ImageProcessor voronoiSkeleton = generateVoronoiSkeletonWithEDM(distanceMap, 4);
        
        // Algorithm 4 - The Bitangent Circle Approach
        // Finally analysis reveals discrete circles are too square-like at small radii.
        //ImageProcessor bitangentSkeleton = generateBitangentSkeletonWithEDM(distanceMap, true);
        // save the skeleton image
        saveFile("skeleton");
        return;
    }

    /**
     * Copies the given image into a new image, and performs the distance
     * map operation on the new image.
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for the original image and must be valid. It is assumed that the given image is binary.
     * <br />Post: On success, a new image of the distance map of the original segmented image is drawn. Also, the local member pixel array is loaded with the original image. On failure, no change is made.  
     * @param ip The Image Processor is required for the original image data.
     * @return The ImageProcessor for the newly created image representing the distance map of the original segmented image. 
     */
    public ImageProcessor generateDistanceMap (ImageProcessor ip) {
        if ((ip == null) || (width < 1) || (height < 1)) {
            return null;
        }
        int x=0, y=0;
        
        String title = "DistanceMap";
        // load the image into memory
        pixel = VascularNetworkToolkit.LoadImage(ip);
        // create a Distance Map of the image       
        ImagePlus dmimp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor dmnp = dmimp.getProcessor();
        for(y=0; y<height; y++) {
            for(x=0; x<width; x++) {
                dmnp.putPixel(x,y,pixel[x][y]); 
            }
        }
        //dmimp.setActivated();
        dmimp.show();
        IJ.selectWindow(title);
        //dmimp.unlock();
        IJ.run("Distance Map");
        //dmimp.lock();
        return dmnp;
    }

    /**
     * Traces the ridges of a given distance map.
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for a distance map and must be valid.
     * <br />Post: On success, a new image of the ridge trace of the distance map of the original segmented image is drawn. Also, the local member pixel array is loaded with the original image. On failure, no change is made.  
     * @param dm The Image Processor is required for the original image data.
     * @return The ImageProcessor for the newly created image representing the ridge trace of the distance map of the original segmented image. 
     * @see #generateDistanceMap
     */
    public ImageProcessor generateRidgeTrace (ImageProcessor dm) {
        if ((dm == null) || (width < 1) || (height < 1)) {
            return null;
        }
        
        String title = "PrethinnedSkeletonization";
        // load the distance map image into memory
        pixel = VascularNetworkToolkit.LoadImage(dm);
        ImagePlus ridgeImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor ridgenp = ridgeImp.getProcessor();

        Coordinate ridgeFinder;
        int x=0, y=0;
        int lightNeighbors, darkNeighbors, matchingNeighbors;
        int lightComponents, fourSquares, maxDelta;
        int lightRegionNeighbors, darkRegionNeighbors, matchingRegionNeighbors;

        /* 
         * Checks the 3 possible cases:
         * 
         * I.   --o  (endpoint)
         * II.  -o-  (ridge/edgepoint)
         * III. >o-  (branch point)
         * IV.  -o-  (branch point
         *      ---   on a 2-thick)
         */
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (pixel[x][y] != WHITE) {
                    ridgeFinder = new Coordinate(x,y,pixel[x][y]);
                    darkNeighbors = dark8Neighbors(ridgeFinder);
                    darkRegionNeighbors = dark24Neighbors(ridgeFinder);
                    lightNeighbors = light8Neighbors(ridgeFinder);
                    matchingNeighbors = matching8Neighbors(ridgeFinder);
                    lightComponents = getLightComponents(ridgeFinder);
                    fourSquares = matching8Neighborhood4Squares(ridgeFinder);
                    lightRegionNeighbors = get5x5Components(LIGHT, ridgeFinder);
/* has lots of noise, not complete
                    if ((lightComponents == 1) && (darkNeighbors <= 4)) {
                        ridgenp.putPixel(x,y,BLACK);
                    }
*/

                    // case 1 -- end points
                    // NOTE: making those 4's threes makes some
                    //       very interesting noise
                    if ((lightComponents == 1) && ((lightNeighbors > 4) || (matchingNeighbors > 4)) ) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                    // ??? - experiment!!
                    else if ((lightComponents == 1) && (matchingNeighbors == 2 && darkNeighbors == 4 && lightNeighbors == 4 && darkRegionNeighbors <= 10)) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                    // ??? - experiment!! newet
                    else if ((lightComponents == 1) && (matchingNeighbors == 3 && darkNeighbors == 5 && lightNeighbors == 3 && darkRegionNeighbors <= 10)) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                    // ULTIMATE EXPERIMENT!1!11!!one11
//                  else if (lightRegionNeighbors >= 2) {
//                      ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
//                  }
                    // case 2 -- edge points
                    else if (lightComponents == 2) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                    // case 3 -- branch points
                    else if (lightComponents >= 3) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                    // case 4 -- thick edge point
                    else if (fourSquares >= 1) {
                        ridgenp.putPixel(x,y,BLACK);
    //                  setColor(x,y,BLACK);
                    }
                }
            }
        }

        // reload the distance map image into memory
        // we kinda broke it ;-)
        //pixel = VascularNetworkToolkit.LoadImage(dm);

        ridgeImp.show();
        IJ.selectWindow(title);
        return ridgenp;     
    }

    /**
     * Performs the Skeletonization (thinning) algorithm on the distance
     * map ridge skeleton.
     * 
     * Sometimes the distance map ridge can be two pixels wide
     * instead of just 1 (imagine a circle with even pixel radius).
     * 
     * Once we have the ridge, thinning will only shift the final solution by
     * 0.5 pixels, so the thinned ridge is still topologically and quantitatively
     * correct, within reasons of error due to the discretization process.
     * 
     * 1) Converts the distance map ridge trace into binary.
     * 2) Performs thinning on the binary image.
     * 3) Uses the thinned binary image as a mask to create a final thinned skeleton.
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for the ridge skeleton and must be valid. In addition, the local member pixel array must be loaded with the distance map.
     * <br />Post: On success, a new image of the thinned distance map ridge trace is drawn. On failure, no change is made.  
     * @param rt The Image Processor is required for the original image data.
     * @return The ImageProcessor for the newly created image representing the thinned distance map ridge trace. 
     * @see #generateDistanceMap
     * @see #generateRidgeTrace
     */
    public ImageProcessor generateThinnedSkeleton (ImageProcessor rt) {
        if ((rt == null) || (width < 1) || (height < 1)) {
            return null;
        }
        int x=0, y=0;
        
        String title = "ThinnedSkeletonization";
        // create new image
        ImagePlus thinImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor thinnp = thinImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                thinnp.putPixel(x,y,rt.getPixel(x,y));
            }
        }

        //thinImp.setActivated();
        thinImp.show();
        IJ.selectWindow(title);
        IJ.run("Skeletonize");

        title = "Skeletonization";
        // create second new window
        ImagePlus skeletonImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor skeletonnp = skeletonImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (thinnp.getPixel(x,y) != WHITE) {
                    skeletonnp.putPixel(x,y,pixel[x][y]);
                }
            }
        }
        skeletonImp.show();
        IJ.selectWindow(title);
        return skeletonnp;      
    }

    /**
     * Performs the thinning algorithm on the 
     * segmented image. The distance map 
     * information is used on the skeleton.
     * 
     * This is a cheap hack way to get a 100% continuous
     * skeleton that does not require any pruning, while
     * preserving many of the advantages of the BMAT
     * approach to the problem.
     * 
     * 1) Converts the distance map into binary.
     * 2) Performs thinning on the binary image.
     * 3) Uses the thinned binary image as a mask 
     *    over the EDM to create a final thinned
     *    skeleton that retains EDM information.
     *    (The skeleton itself is NOT a true BMAT.)
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for the ridge skeleton and must be valid. In addition, the local member pixel array must be loaded with the distance map.
     * <br />Post: On success, a new image of the thinned distance map ridge trace is drawn. On failure, no change is made.  
     * @param rt The Image Processor is required for the original image data.
     * @return The ImageProcessor for the newly created image representing the thinned distance map ridge trace. 
     * @see #generateDistanceMap
     */
    public ImageProcessor generateThinningSkeletonWithEDM (ImageProcessor dm) {
        if (dm == null) {
            return null;
        }
        int x=0, y=0;
        
        // load the distance map image into memory
        pixel = VascularNetworkToolkit.LoadImage(dm);   
        
        String title = "Skeletonization";
        // create new image
        ImagePlus thinImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor thinnp = thinImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (getColor(x,y) != WHITE) {
                    thinnp.putPixel(x,y,BLACK);
                }
            }
        }

        //thinImp.setActivated();
        thinImp.show();
        IJ.selectWindow(title);
        IJ.run("Skeletonize");
        IJ.wait(100);

        title = "EDMSkeletonization";
        // create second new window
        ImagePlus skeletonImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor skeletonnp = skeletonImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (thinnp.getPixel(x,y) != WHITE) {
                    skeletonnp.putPixel(x,y,getColor(x,y));
                }
            }
        }
        skeletonImp.show();
        IJ.selectWindow(title);
        IJ.wait(100);
        return skeletonnp;      
    }
    
    /**
     * Performs the fancy pants Voronoi/Delaunay algorithm stuff
     * that Paul Chew, chew@cs.cornell.edu wrote.
     * I got his code from http://www.cs.cornell.edu/Info/People/chew/Delaunay.html
     *  
     * 1) Finds the edges of the distance map and inputs them into the Voronoi algorithm.
     * 2) Draws the Voronoi region boundaries.
     * 3) Uses the boundaries as a mask 
     *    over the EDM to create a final thinned
     *    skeleton that retains EDM information.
     *    (The skeleton itself may NOT be a true BMAT.)
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for the ridge skeleton and must be valid. In addition, the local member pixel array must be loaded with the distance map.
     * <br />Post: On success, a new image of the thinned distance map ridge trace is drawn. On failure, no change is made.  
     * @param rt The Image Processor is required for the original image data.
     * @return The ImageProcessor for the newly created image representing the thinned distance map ridge trace. 
     * @see #generateDistanceMap
     */
    public ImageProcessor generateVoronoiSkeletonWithEDM (ImageProcessor dm, int searchMax) {
        if (dm == null) {
            return null;
        }

        int x=0, y=0;
        String title;
        
        // load the distance map image into memory
        pixel = VascularNetworkToolkit.LoadImage(dm);   
        
        title = "Skeletonization";
        // create new image
        ImagePlus thinImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor thinnp = thinImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (getColor(x,y) != WHITE) {
                    thinnp.putPixel(x,y,BLACK);
                }
            }
        }
        //thinImp.setActivated();
        thinImp.show();
        IJ.selectWindow(title);
        IJ.run("Skeletonize");

        DelaunayTriangulation dt;   // The Delaunay triangulation
        Simplex initialTriangle;        // The large initial triangle
        Pnt point;                          // the point datastructure
        int initialSize = 10000;      // Controls size of initial triangle
        int searched, deltaX, deltaY, segmentPoints;
        
        // prepare the Voronoi/Delaunay algorithm!!
        initialTriangle = new Simplex(new Pnt[] {
            new Pnt(-initialSize, -initialSize),
            new Pnt( initialSize, -initialSize),
            new Pnt(           0,  initialSize)});
        // search stack to travel the pixels
        Stack searchStack = new Stack();
        Coordinate seeker;
        boolean found;
        
        title = "VoronoiDelaunay";
        // create new image
        ImagePlus voronoiImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor voronoinp = voronoiImp.getProcessor();

        // load the distance map image into memory
        pixel = VascularNetworkToolkit.LoadImage(dm);

        // clear the datastructure
        dt = new DelaunayTriangulation(initialTriangle);
        // use the distance map to fill the delaunay structure  
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (getColor(x,y) == WHITE-1) { // if at border

                    searchStack.push(new Coordinate(x,y,WALL));
                    searched = 0;
                    segmentPoints = 0;
                    while (!searchStack.isEmpty()) {
                        seeker = (Coordinate)searchStack.pop();
                        if (getColor(seeker.getX(),seeker.getY()) != WHITE) {
                            if (searched <= 0) {
                                // place a point
                                segmentPoints ++;
                                point = new Pnt(seeker.getX(), seeker.getY());
                                dt.delaunayPlace(point);
                                searched = searchMax;                               
                            }
                            searched --;
                            setColor(seeker.getX(),seeker.getY(),WHITE);
                            // this is terrible but it ensures only one pixel (if any) are added
                            found = false;
                            deltaX=1;deltaY=0; // to the right
                            if (getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }
                            deltaX=0;deltaY=-1; // up
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=-1;deltaY=0; // to the left
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=0;deltaY=1; // down
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=1;deltaY=-1; // upper right
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=-1;deltaY=-1; // upper left
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=-1;deltaY=1; // lower left
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                            deltaX=1;deltaY=1; // lower right
                            if (!found && getColor(seeker.getX()+deltaX,seeker.getY()+deltaY) == WHITE-1) {
                                searchStack.push(new Coordinate(seeker.getX()+deltaX,seeker.getY()+deltaY,WALL));
                                found = true;
                            }  
                        }
                    }
                    System.out.println("There were "+segmentPoints+" points found for a shape at "+x+","+y);
                    // if it was a too-small-blob, use the thinning algorithm
                    // draw it (either voronoi or thinning)
                    if (segmentPoints <= 1 && false ) { // skip this for now
                        searchStack.push(new Coordinate(x,y,WALL));
                        while (!searchStack.isEmpty()) {
                            seeker = (Coordinate)searchStack.pop();
                            //if 
                        }
                    } else {
                        // draw code should go here
                    }
                }
            }
        }
        
        // draw all Voronoi boundaries (including outside ones we don't want)
        // Loop through all the edges of the DT (each is done twice)
        Simplex triangle, other;
        Pnt p,q;
        for (Iterator it = dt.iterator(); it.hasNext();) {
            triangle = (Simplex) it.next();
            for (Iterator otherIt = dt.neighbors(triangle).iterator(); otherIt.hasNext();) {
                other = (Simplex) otherIt.next();
                p = Pnt.circumcenter((Pnt[]) triangle.toArray(new Pnt[0]));
                q = Pnt.circumcenter((Pnt[]) other.toArray(new Pnt[0]));
                drawLine(voronoinp, (int)p.coord(0), (int)p.coord(1), (int)q.coord(0), (int)q.coord(1), BLACK);
            }
        }   
        /*for (Iterator it = dt.iterator(); it.hasNext();) {
            for (Iterator otherIt = ((Simplex) it.next()).iterator(); otherIt.hasNext();) {
                p = (Pnt) otherIt.next();
                drawCircle(voronoinp, (int)p.coord(0),(int)p.coord(1), 2, BLACK);
            }
        }*/                             
    
        
        // reload the distance map image into memory
        pixel = VascularNetworkToolkit.LoadImage(dm);

        // now use memory to delete any lines drawn on the background
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (getColor(x,y) == WHITE && voronoinp.getPixel(x,y) != WHITE) {
                    voronoinp.putPixel(x,y,WHITE);
                }
            }
        }

        voronoiImp.show();
        IJ.selectWindow(title);
        //IJ.run("Skeletonize");

        title = "EDMVoronoiSkeletonization";
        // create second new window
        ImagePlus skeletonImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor skeletonnp = skeletonImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (voronoinp.getPixel(x,y) != WHITE) {
                    skeletonnp.putPixel(x,y,getColor(x,y));
                }
            }
        }
        skeletonImp.show();
        IJ.selectWindow(title);
        return skeletonnp;      
    }
    
    /**
     * Calculates the distance between two radian angles.
     * Only slightly trickier than it sounds!
     * 
     * <p>Pre: None.
     * <br />Post: None. Static method, after all..
     * @param angleA The first angle. Must be in radians.
     * @param angleB The second angle. Must be in radians.
     * @return Returns the radian distance between the two angles. 
     */
    public static double calculateAngleDistance(double angleA, double angleB) {
        double finalAngle = angleA - angleB;
        // make sure our given angles are within bounds 0-360
        while (angleA < 0) {
            angleA += 2*Math.PI; // add 360 degrees until not negative
        }
        while (angleB < 0) {
            angleB += 2*Math.PI; // add 360 degrees until not negative
        }
        while (angleA >= 2*Math.PI) {
            angleA -= 2*Math.PI; // subtract 360 degrees until less than 360
        }
        while (angleB >= 2*Math.PI) {
            angleB -= 2*Math.PI; // subtract 360 degrees until less than 360
        }
        // compute our angle difference
        finalAngle = angleA - angleB;
        // make sure our angle difference is within bounds
        while (finalAngle < 0) {
            finalAngle += 2*Math.PI; // add 360 degrees until not negative
        }
        while (finalAngle >= 2*Math.PI) {
            finalAngle -= 2*Math.PI; // subtract 360 degrees until lessthan 360
        }
        // if angle is greater than 180, then take 360-angle
        if (finalAngle > Math.PI) { // not a typo
            finalAngle = 2*Math.PI - finalAngle;
        }
        // our angle is now between 0 and 180, a good difference
        return finalAngle;
    }

    /**
     * Performs a search for the BMAT by testing for 
     * maximal radius bitangent circles that have no 
     * part outside of the object. The envelope of
     * the centers of these circles is the ridge.
     * 
     * <p>Pre: The dimensions width and height must have been loaded prior to calling this method. The given ImageProcessor must be for the EDM and must be valid.
     * <br />Post: On success, a new image of the thinned distance map ridge trace is drawn. On failure, no change is made.  
     * @param dm The Image Processor is required for the EDM image data.
     * @return The ImageProcessor for the newly created image representing the thinned distance map ridge trace. 
     * @see #generateDistanceMap
     */
    public ImageProcessor generateBitangentSkeletonWithEDM (ImageProcessor dm) {
        if (dm == null) {
            return null;
        }

        int x=0, y=0,dx=0,dy=0,lastx=0,lasty=0;
        int radius=0,borderPixels=0;
        
        double angle=0;
        double maxAngle = Math.toRadians(60);
        String title;
        
        title = "Scratch Paper";
        // create new image
        ImagePlus scratchImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor scratchnp = scratchImp.getProcessor();
        
        if (animatedDisplay) {
            for (y=0;y<height;y++) {
                for (x=0;x<width;x++) {
                    if (dm.getPixel(x,y) != WHITE) {
                        scratchnp.putPixel(x,y,dm.getPixel(x,y)-100);
                    }
                }
            }
        }
        
        
        title = "Envelope of Centers of Bitangent Circles";
        // create new image
        ImagePlus bitangentImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor bitangentnp = bitangentImp.getProcessor();

        if (animatedDisplay) {          
            for (y=0; y<height; y++) {
                for (x=0; x<width; x++) {
                    if (dm.getPixel(x,y) != WHITE) {
                        bitangentnp.putPixel(x,y,BLACK);
                    }
                }
            }
            
            scratchImp.show();
            bitangentImp.show();
        }
        
        // ??? actually compute everything
        for (y=0; y<height; y++) {
            IJ.showProgress(VascularNetworkToolkit.convertPercentage(y,height)); 
            for (x=0; x<width; x++) {
                // look for non-white pixels on the distance map
                radius = WHITE - dm.getPixel(x,y);
                // draw a circle and see if it's a maximal bitangent
                // circle without overlap
                if (radius != BLACK) {                  
                    // draw a circle
                    VascularNetworkToolkit.drawCircle(scratchnp,x,y,radius-2,BLACK);
                    if (animatedDisplay) {
                        for (dy=y-radius*2;dy<=y+radius*2;dy++) {
                            for (dx=x-radius*2;dx<=x+radius*2;dx++) {
                                if (dx >= 0 && dy>= 0 && dx < width && dy < height) {
                                    if (scratchnp.getPixel(dx,dy) != BLACK && dm.getPixel(dx,dy) != WHITE) {
                                        scratchnp.putPixel(dx,dy,dm.getPixel(dx,dy)-100);
                                    }
                                }
                            }                   
                        }
                    }
                    // count the results
                    borderPixels = 0;
                    for (dy=y-radius;dy<=y+radius;dy++) {
                        for (dx=x-radius;dx<=x+radius;dx++) {
                            if (scratchnp.getPixel(dx,dy) == BLACK) {
                                if (dm.getPixel(dx,dy) <= WHITE-1 && dm.getPixel(dx,dy) >= WHITE-5) {
                                    if (borderPixels == 0) {
                                        // tangent at this point
                                        borderPixels ++;
                                        lastx = dx;
                                        lasty = dy;                                 
                                        angle = Math.atan2((double)(dy-y),(double)(dx-x));
                                        if (animatedDisplay) {
                                            VascularNetworkToolkit.drawLine(scratchnp,dx,dy,x,y,WHITE);
                                        }
                                    } else if (borderPixels == 1 && (Math.abs(dx-lastx)+Math.abs(dy-lasty) > Math.sqrt((double)(2*radius))) && (calculateAngleDistance(angle, Math.atan2((double)(dy-y),(double)(dx-x))) > maxAngle)) {
                                        borderPixels ++;
                                        if (animatedDisplay) {
                                            VascularNetworkToolkit.drawLine(scratchnp,dx,dy,x,y,WHITE);
                                        }
                                    }
                                } else if (dm.getPixel(dx,dy) == WHITE) {
                                    // overlap, failure. force abort
                                    dy = y + radius;
                                    dx = x + radius;
                                    borderPixels = 0;
                                }
                            }
                        }
                    }

                    // conclusion: only draw ridge points
                    if (borderPixels >= 2) { // at least bitangent! success!
                        bitangentnp.putPixel(x,y,BLACK);
                        if (animatedDisplay) {
                            scratchImp.updateAndDraw();
                            bitangentImp.updateAndDraw();
                        }
                    } else if (animatedDisplay) {
                        bitangentnp.putPixel(x,y,WHITE);                        
                        bitangentImp.updateAndDraw();
                    }
                    // erase the circle
                    if (animatedDisplay) {
                        for (dy=y-radius*2;dy<=y+radius*2;dy++) {
                            for (dx=x-radius*2;dx<=x+radius*2;dx++) {
                                if (dx >= 0 && dy>= 0 && dx < width && dy < height) {
                                    if (dm.getPixel(dx,dy) == WHITE) {
                                        scratchnp.putPixel(dx,dy,WHITE);                                        
                                    } else {
                                        scratchnp.putPixel(dx,dy,dm.getPixel(dx,dy)-100);
                                    }
                                }
                            }                   
                        }
                    } else {
                        VascularNetworkToolkit.drawCircle(scratchnp,x,y,radius,WHITE);                      
                    }
                    
                }
            }           
        }

        // reveal the wonders to the world!
        bitangentImp.show();
        IJ.selectWindow(title);
        //IJ.run("Skeletonize");

        title = "EDMBitangentSkeletonization";
        // create second new window
        ImagePlus skeletonImp = NewImage.createByteImage (title, width, height, 1, NewImage.FILL_WHITE);
        ImageProcessor skeletonnp = skeletonImp.getProcessor();
        for (x=0; x<width; x++) {
            for (y=0; y<height; y++) {
                if (bitangentnp.getPixel(x,y) != WHITE) {
                    skeletonnp.putPixel(x,y,dm.getPixel(x,y));
                }
            }
        }
        skeletonImp.show();
        IJ.selectWindow(title);
        return skeletonnp;      
    }
}