From fiquère, 9 Months ago, written in Plain Text.
This paste will die in 1 Second.
Embed
  1. package Hough_Package;
  2.  
  3. /** Hough_Circle.java:
  4.  
  5.  This program is free software; you can redistribute it and/or modify
  6.  it under the terms of the GNU General Public License as published by
  7.  the Free Software Foundation; either version 2 of the License, or
  8.  (at your option) any later version.
  9.  
  10.  This program is distributed in the hope that it will be useful,
  11.  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13.  GNU General Public License for more details.
  14.  
  15.  @author Ben Smith (benjamin.smith@berkeley.edu)
  16.  @Based on original plugin implementation by Hemerson Pistori (pistori@ec.ucdb.br) and Eduardo Rocha Costa
  17.  @created February 4, 2017
  18.  
  19.  The Hough Transform implementation was based on
  20.  Mark A. Schulze applet (http://www.markschulze.net/)
  21.  
  22. */
  23.  
  24. import ij.*;
  25. import ij.measure.Calibration;
  26. import ij.process.*;
  27. import java.awt.*;
  28. import ij.plugin.HyperStackConverter;
  29. import ij.plugin.frame.*;
  30. import java.util.concurrent.atomic.AtomicInteger;
  31. import javax.swing.*;
  32. import ij.measure.ResultsTable;
  33. import ij.plugin.filter.Analyzer;
  34. import java.util.Arrays;
  35. import java.util.List;
  36.  
  37. /**
  38.  *
  39.  * @author Ben
  40.  */
  41. public class Hough_Circle extends SwingWorker<Integer, String>{
  42.     //Input parameters
  43.     private int radiusMin;  // Find circles with radius grater or equal radiusMin - argument syntax: "min=#"
  44.     private int radiusMax;  // Find circles with radius less or equal radiusMax - argument syntax: "max=#"
  45.     private int radiusInc;  // Increment used to go from radiusMin to radiusMax - argument syntax: "inc=#"
  46.     private int minCircles; // Minumum number of circles to be found - argument syntax: "minCircles=#"    
  47.     private int maxCircles; // Maximum number of circles to be found - argument syntax: "maxCircles=#"
  48.     private int threshold = -1; // An alternative to maxCircles. All circles with a value in the hough space greater then threshold are marked. Higher thresholds result in fewer circles. - argument syntax: "threshold=#"
  49.     private double thresholdRatio; //Ratio input from GUI that expresses threshold as ratio of resolution (highest possible # of votes)
  50.     private int resolution; //The number of steps to use per transform (i.e. number of voting rounds)
  51.     private double ratio; // Ratio of found circle radius to clear out surrounding neighbors
  52.     private int searchBand = 0; //The +/- range of radii to search for relative to the last found radius - argument syntax: "bandwidth=#"
  53.     private int searchRadius = 0; //The search radius to look for the next centroid relative to the last found centroid - argument syntax: "radius=#"
  54.     private boolean reduce = false; //Cap the transform resolution by removeing redundant steps
  55.     private boolean local = false; //Whether or not the search is going to be local
  56.    
  57.     //Output parameters
  58.     private boolean houghSeries = false; //Contains whether the user wants the Hough series stack as an output - argument syntax: "show_raw"
  59.     private boolean showCircles = false; //Contains whether the user wants the circles found as an output - argument syntax: "show_mask"
  60.     private boolean showID = false; //Contains whether the user wants a map of centroids and radii outputed from search - argument syntax: "show_centroids"
  61.     private boolean showScores = false; //Contains whether the user wants a map of centroids and Hough scores outputed from search - argument syntax: "show_scores"
  62.     private boolean results = false; //Contains whether the user wants to export the measurements to a reuslts table
  63.     private Calibration pixelCal;
  64.     private String pixelUnits; //Stores the unit of measurement fo the pixels
  65.     private double pixelDimensions; //Stores the size of each pixel
  66.     private int[] imageDimensions; //Pixel dimensions of stack
  67.     private String timeUnits; //Frame units
  68.     private double timeDimension; //Time step per frame
  69.    
  70.     private String currentStatus = ""; //String for outputting current status
  71.     private boolean isGUI; //Whether a GUI is active (or a macro called the plugin)
  72.     private final static int GUI_UPDATE_DELAY = 100; //How long to wait between GUI updates
  73.     private boolean cancelThread = false;
  74.            
  75.     //Hough transform variables
  76.     private ImagePlus imp; //Initalize the variable to hold the image
  77.     private boolean isStack = false; //True if there is more than one slice in the input data
  78.     private int stackSlices; //number of slices in the stack
  79.     private int maxHough; //Contains the brights pixel in the entire Hough array
  80.     private Point maxPoint; //Stores the location of the brightest pixel in the Hough array
  81.     private int maxRadius; //Stores the radius of the brightest pixel in the Hough array
  82.     private Rectangle r; //Stores the ROI on the original image
  83.     private float imageValues[]; // Raw image (returned by ip.getPixels()) - float is used to allow 8, 16 or 32 bit images
  84.     private int houghValues[][][]; // Hough Space Values [X coord][Y coord][radius index]
  85.     private int localHoughValues[][][][]; //Local Hough space [circle#][X coord][Y coord][radius index]
  86.     private int localHoughParameters[][]; //Array to pass local Hough space parameters to centroid search [circle#][parameter vector]
  87.     private int width; // ROI width
  88.     private int height;  // ROI height
  89.     private int depth;  // Number of slices
  90.     private int fullWidth; // Image Width
  91.     private int fullHeight; //Image Height
  92.     private int offx;   // ROI x origin position
  93.     private int offy;   // ROI y origin position
  94.     private Point centerPoint[]; // Center Points of the Circles Found.
  95.     private int centerRadii[]; //Corresponding radii of the cricles marked by the center points
  96.     private int houghScores[]; //Corresponding Hough scores for each centroid
  97.     private int circleID[]; //Corresponding ID # for each centroid
  98.     private int circleIDcounter; //Counter for keeping track of current max ID #
  99.     private int lut[][][]; // LookUp Table for x and y tranform shifts in an octahedral manner
  100.     private int lutSize; //Stores the actual number of transforms performed (<=selected resolution)
  101.     private int nCircles; //Number of circles found during search - <= maxCircles
  102.     private int nCirlcesPrev; //Stores nCirlces from last iteration
  103.     private boolean localSearch = false; //Record whether a local-only search was done for the frame
  104.    
  105.     //Variables for storing the results and exporting result images
  106.     private ImagePlus houghPlus;
  107.     private ImageStack houghStack;
  108.     private ImagePlus circlePlus;
  109.     private ImageStack circleStack;
  110.     private ImagePlus idPlus;
  111.     private ImageStack idStack;
  112.     private ImagePlus scorePlus;
  113.     private ImageStack scoreStack;
  114.     private ImageProcessor circlesip;
  115.     private ImageProcessor IDip;
  116.     private ImageProcessor scoresip;
  117.     private ResultsTable rt;
  118.     private String method;
  119.    
  120.     //Build LUTs for colorimetric output
  121.     private final static byte[] RAND_R = new byte[] {(byte) 0, (byte) 231, (byte) 25, (byte) 27, (byte) 112, (byte) 89, (byte) 53, (byte) 77, (byte) 155, (byte) 153, (byte) 91, (byte) 19, (byte) 25, (byte) 135, (byte) 187, (byte) 40, (byte) 52, (byte) 39, (byte) 147, (byte) 72, (byte) 92, (byte) 23, (byte) 5, (byte) 202, (byte) 52, (byte) 50, (byte) 108, (byte) 66, (byte) 238, (byte) 46, (byte) 136, (byte) 8, (byte) 162, (byte) 126, (byte) 75, (byte) 60, (byte) 135, (byte) 167, (byte) 136, (byte) 244, (byte) 161, (byte) 204, (byte) 158, (byte) 112, (byte) 192, (byte) 2, (byte) 98, (byte) 207, (byte) 233, (byte) 121, (byte) 191, (byte) 18, (byte) 210, (byte) 151, (byte) 252, (byte) 155, (byte) 139, (byte) 42, (byte) 188, (byte) 170, (byte) 127, (byte) 151, (byte) 38, (byte) 16, (byte) 227, (byte) 120, (byte) 137, (byte) 2, (byte) 72, (byte) 71, (byte) 214, (byte) 130, (byte) 194, (byte) 96, (byte) 102, (byte) 33, (byte) 92, (byte) 47, (byte) 211, (byte) 226, (byte) 237, (byte) 132, (byte) 244, (byte) 45, (byte) 82, (byte) 73, (byte) 239, (byte) 198, (byte) 118, (byte) 191, (byte) 170, (byte) 238, (byte) 177, (byte) 246, (byte) 66, (byte) 5, (byte) 44, (byte) 3, (byte) 18, (byte) 40, (byte) 226, (byte) 94, (byte) 130, (byte) 216, (byte) 68, (byte) 89, (byte) 150, (byte) 40, (byte) 96, (byte) 155, (byte) 95, (byte) 5, (byte) 132, (byte) 88, (byte) 124, (byte) 241, (byte) 143, (byte) 15, (byte) 1, (byte) 220, (byte) 46, (byte) 220, (byte) 216, (byte) 155, (byte) 190, (byte) 4, (byte) 98, (byte) 33, (byte) 142, (byte) 153, (byte) 185, (byte) 61, (byte) 153, (byte) 160, (byte) 164, (byte) 102, (byte) 75, (byte) 8, (byte) 26, (byte) 231, (byte) 160, (byte) 245, (byte) 119, (byte) 69, (byte) 95, (byte) 137, (byte) 232, (byte) 249, (byte) 89, (byte) 220, (byte) 54, (byte) 3, (byte) 43, (byte) 134, (byte) 172, (byte) 31, (byte) 106, (byte) 234, (byte) 132, (byte) 1, (byte) 29, (byte) 65, (byte) 12, (byte) 107, (byte) 52, (byte) 196, (byte) 95, (byte) 165, (byte) 4, (byte) 101, (byte) 124, (byte) 240, (byte) 203, (byte) 50, (byte) 142, (byte) 231, (byte) 5, (byte) 23, (byte) 43, (byte) 198, (byte) 41, (byte) 80, (byte) 66, (byte) 232, (byte) 86, (byte) 151, (byte) 126, (byte) 51, (byte) 122, (byte) 111, (byte) 162, (byte) 192, (byte) 148, (byte) 210, (byte) 5, (byte) 251, (byte) 40, (byte) 176, (byte) 32, (byte) 2, (byte) 22, (byte) 91, (byte) 180, (byte) 31, (byte) 218, (byte) 9, (byte) 10, (byte) 178, (byte) 27, (byte) 115, (byte) 54, (byte) 199, (byte) 163, (byte) 10, (byte) 180, (byte) 222, (byte) 102, (byte) 198, (byte) 202, (byte) 220, (byte) 96, (byte) 145, (byte) 237, (byte) 127, (byte) 172, (byte) 51, (byte) 64, (byte) 47, (byte) 84, (byte) 37, (byte) 196, (byte) 147, (byte) 244, (byte) 206, (byte) 244, (byte) 43, (byte) 170, (byte) 186, (byte) 162, (byte) 157, (byte) 125, (byte) 230, (byte) 212, (byte) 140, (byte) 186, (byte) 127, (byte) 12, (byte) 48, (byte) 0, (byte) 92, (byte) 168, (byte) 73, (byte) 228, (byte) 121, (byte) 248, (byte) 255};
  122.     private final static byte[] RAND_G = new byte[] {(byte) 0, (byte) 217, (byte) 131, (byte) 18, (byte) 66, (byte) 198, (byte) 59, (byte) 170, (byte) 108, (byte) 8, (byte) 167, (byte) 207, (byte) 248, (byte) 138, (byte) 215, (byte) 204, (byte) 95, (byte) 249, (byte) 102, (byte) 66, (byte) 1, (byte) 138, (byte) 71, (byte) 33, (byte) 209, (byte) 80, (byte) 72, (byte) 87, (byte) 129, (byte) 24, (byte) 200, (byte) 98, (byte) 16, (byte) 160, (byte) 47, (byte) 1, (byte) 224, (byte) 49, (byte) 126, (byte) 41, (byte) 235, (byte) 205, (byte) 84, (byte) 249, (byte) 222, (byte) 125, (byte) 5, (byte) 58, (byte) 90, (byte) 137, (byte) 73, (byte) 98, (byte) 109, (byte) 167, (byte) 91, (byte) 106, (byte) 68, (byte) 23, (byte) 10, (byte) 241, (byte) 51, (byte) 35, (byte) 77, (byte) 112, (byte) 243, (byte) 235, (byte) 53, (byte) 41, (byte) 21, (byte) 134, (byte) 106, (byte) 60, (byte) 231, (byte) 107, (byte) 124, (byte) 96, (byte) 131, (byte) 62, (byte) 213, (byte) 220, (byte) 78, (byte) 96, (byte) 44, (byte) 106, (byte) 128, (byte) 96, (byte) 129, (byte) 164, (byte) 59, (byte) 86, (byte) 45, (byte) 87, (byte) 237, (byte) 192, (byte) 171, (byte) 93, (byte) 223, (byte) 206, (byte) 158, (byte) 42, (byte) 255, (byte) 175, (byte) 142, (byte) 177, (byte) 48, (byte) 67, (byte) 221, (byte) 123, (byte) 125, (byte) 182, (byte) 135, (byte) 31, (byte) 163, (byte) 211, (byte) 220, (byte) 28, (byte) 81, (byte) 113, (byte) 235, (byte) 236, (byte) 30, (byte) 104, (byte) 238, (byte) 227, (byte) 18, (byte) 134, (byte) 193, (byte) 15, (byte) 122, (byte) 151, (byte) 6, (byte) 132, (byte) 0, (byte) 128, (byte) 88, (byte) 20, (byte) 130, (byte) 87, (byte) 127, (byte) 55, (byte) 13, (byte) 76, (byte) 15, (byte) 124, (byte) 251, (byte) 160, (byte) 177, (byte) 179, (byte) 14, (byte) 7, (byte) 189, (byte) 49, (byte) 210, (byte) 77, (byte) 63, (byte) 132, (byte) 17, (byte) 84, (byte) 214, (byte) 63, (byte) 206, (byte) 252, (byte) 184, (byte) 168, (byte) 109, (byte) 129, (byte) 89, (byte) 70, (byte) 66, (byte) 151, (byte) 151, (byte) 100, (byte) 27, (byte) 123, (byte) 35, (byte) 210, (byte) 112, (byte) 183, (byte) 135, (byte) 97, (byte) 19, (byte) 7, (byte) 159, (byte) 221, (byte) 46, (byte) 253, (byte) 42, (byte) 229, (byte) 101, (byte) 187, (byte) 10, (byte) 33, (byte) 152, (byte) 138, (byte) 236, (byte) 208, (byte) 117, (byte) 230, (byte) 98, (byte) 200, (byte) 139, (byte) 200, (byte) 255, (byte) 19, (byte) 219, (byte) 183, (byte) 181, (byte) 3, (byte) 43, (byte) 143, (byte) 130, (byte) 75, (byte) 21, (byte) 228, (byte) 121, (byte) 208, (byte) 51, (byte) 82, (byte) 41, (byte) 233, (byte) 56, (byte) 220, (byte) 190, (byte) 176, (byte) 136, (byte) 108, (byte) 88, (byte) 118, (byte) 228, (byte) 206, (byte) 11, (byte) 170, (byte) 236, (byte) 232, (byte) 204, (byte) 241, (byte) 241, (byte) 69, (byte) 71, (byte) 28, (byte) 143, (byte) 207, (byte) 52, (byte) 188, (byte) 183, (byte) 80, (byte) 4, (byte) 222, (byte) 162, (byte) 30, (byte) 213, (byte) 228, (byte) 119, (byte) 142, (byte) 10, (byte) 255};
  123.     private final static byte[] RAND_B = new byte[] {(byte) 0, (byte) 43, (byte) 203, (byte) 219, (byte) 253, (byte) 158, (byte) 5, (byte) 190, (byte) 72, (byte) 11, (byte) 23, (byte) 233, (byte) 220, (byte) 246, (byte) 53, (byte) 10, (byte) 90, (byte) 49, (byte) 215, (byte) 182, (byte) 74, (byte) 50, (byte) 27, (byte) 107, (byte) 39, (byte) 48, (byte) 192, (byte) 134, (byte) 247, (byte) 89, (byte) 14, (byte) 3, (byte) 67, (byte) 65, (byte) 30, (byte) 136, (byte) 78, (byte) 129, (byte) 178, (byte) 138, (byte) 186, (byte) 204, (byte) 5, (byte) 160, (byte) 18, (byte) 103, (byte) 255, (byte) 162, (byte) 42, (byte) 128, (byte) 213, (byte) 204, (byte) 49, (byte) 80, (byte) 181, (byte) 130, (byte) 60, (byte) 185, (byte) 31, (byte) 203, (byte) 184, (byte) 89, (byte) 108, (byte) 190, (byte) 109, (byte) 157, (byte) 231, (byte) 62, (byte) 128, (byte) 96, (byte) 150, (byte) 153, (byte) 160, (byte) 54, (byte) 24, (byte) 143, (byte) 90, (byte) 216, (byte) 128, (byte) 120, (byte) 87, (byte) 244, (byte) 15, (byte) 213, (byte) 235, (byte) 142, (byte) 140, (byte) 33, (byte) 98, (byte) 164, (byte) 202, (byte) 38, (byte) 84, (byte) 63, (byte) 229, (byte) 163, (byte) 28, (byte) 239, (byte) 210, (byte) 131, (byte) 79, (byte) 83, (byte) 168, (byte) 79, (byte) 89, (byte) 170, (byte) 26, (byte) 168, (byte) 149, (byte) 217, (byte) 31, (byte) 20, (byte) 11, (byte) 141, (byte) 121, (byte) 139, (byte) 21, (byte) 58, (byte) 194, (byte) 75, (byte) 110, (byte) 234, (byte) 16, (byte) 126, (byte) 24, (byte) 41, (byte) 62, (byte) 88, (byte) 232, (byte) 38, (byte) 243, (byte) 83, (byte) 195, (byte) 58, (byte) 84, (byte) 106, (byte) 151, (byte) 32, (byte) 146, (byte) 24, (byte) 140, (byte) 217, (byte) 176, (byte) 186, (byte) 174, (byte) 170, (byte) 250, (byte) 1, (byte) 48, (byte) 141, (byte) 99, (byte) 213, (byte) 127, (byte) 46, (byte) 97, (byte) 12, (byte) 143, (byte) 237, (byte) 153, (byte) 185, (byte) 72, (byte) 170, (byte) 23, (byte) 222, (byte) 216, (byte) 23, (byte) 92, (byte) 232, (byte) 54, (byte) 64, (byte) 72, (byte) 128, (byte) 182, (byte) 38, (byte) 192, (byte) 138, (byte) 115, (byte) 162, (byte) 231, (byte) 2, (byte) 143, (byte) 117, (byte) 167, (byte) 196, (byte) 4, (byte) 182, (byte) 220, (byte) 52, (byte) 252, (byte) 34, (byte) 2, (byte) 233, (byte) 132, (byte) 135, (byte) 230, (byte) 36, (byte) 199, (byte) 228, (byte) 222, (byte) 43, (byte) 119, (byte) 7, (byte) 148, (byte) 59, (byte) 10, (byte) 35, (byte) 158, (byte) 237, (byte) 17, (byte) 116, (byte) 110, (byte) 129, (byte) 172, (byte) 233, (byte) 172, (byte) 47, (byte) 114, (byte) 26, (byte) 115, (byte) 212, (byte) 177, (byte) 157, (byte) 74, (byte) 183, (byte) 102, (byte) 132, (byte) 151, (byte) 18, (byte) 242, (byte) 242, (byte) 12, (byte) 138, (byte) 21, (byte) 159, (byte) 165, (byte) 213, (byte) 230, (byte) 147, (byte) 174, (byte) 206, (byte) 116, (byte) 9, (byte) 242, (byte) 233, (byte) 202, (byte) 205, (byte) 116, (byte) 169, (byte) 80, (byte) 53, (byte) 161, (byte) 239, (byte) 211, (byte) 73, (byte) 96, (byte) 255};
  124.     private final static LUT RAND_LUT = new LUT(RAND_R, RAND_G, RAND_B);
  125.    
  126.     private final static byte[] GYR_R = new byte[] {(byte) 0, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 253, (byte) 251, (byte) 249, (byte) 247, (byte) 245, (byte) 243, (byte) 241, (byte) 239, (byte) 237, (byte) 235, (byte) 233, (byte) 231, (byte) 229, (byte) 227, (byte) 225, (byte) 223, (byte) 221, (byte) 219, (byte) 217, (byte) 215, (byte) 213, (byte) 211, (byte) 209, (byte) 207, (byte) 205, (byte) 203, (byte) 201, (byte) 199, (byte) 197, (byte) 195, (byte) 193, (byte) 191, (byte) 189, (byte) 187, (byte) 185, (byte) 183, (byte) 181, (byte) 179, (byte) 177, (byte) 175, (byte) 173, (byte) 171, (byte) 169, (byte) 167, (byte) 165, (byte) 163, (byte) 161, (byte) 159, (byte) 157, (byte) 155, (byte) 153, (byte) 151, (byte) 149, (byte) 147, (byte) 145, (byte) 143, (byte) 141, (byte) 139, (byte) 137, (byte) 135, (byte) 133, (byte) 131, (byte) 129, (byte) 127, (byte) 125, (byte) 123, (byte) 121, (byte) 119, (byte) 117, (byte) 115, (byte) 113, (byte) 111, (byte) 109, (byte) 107, (byte) 105, (byte) 103, (byte) 101, (byte) 99, (byte) 97, (byte) 95, (byte) 93, (byte) 91, (byte) 89, (byte) 87, (byte) 85, (byte) 83, (byte) 81, (byte) 79, (byte) 77, (byte) 75, (byte) 73, (byte) 71, (byte) 69, (byte) 67, (byte) 65, (byte) 63, (byte) 61, (byte) 59, (byte) 57, (byte) 55, (byte) 53, (byte) 51, (byte) 49, (byte) 47, (byte) 45, (byte) 43, (byte) 41, (byte) 39, (byte) 37, (byte) 35, (byte) 33, (byte) 31, (byte) 29, (byte) 27, (byte) 25, (byte) 23, (byte) 21, (byte) 19, (byte) 17, (byte) 15, (byte) 13, (byte) 11, (byte) 9, (byte) 7, (byte) 5, (byte) 3, (byte) 1};
  127.     private final static byte[] GYR_G = new byte[] {(byte) 0, (byte) 1, (byte) 3, (byte) 5, (byte) 7, (byte) 9, (byte) 11, (byte) 13, (byte) 15, (byte) 17, (byte) 19, (byte) 21, (byte) 23, (byte) 25, (byte) 27, (byte) 29, (byte) 31, (byte) 33, (byte) 35, (byte) 37, (byte) 39, (byte) 41, (byte) 43, (byte) 45, (byte) 47, (byte) 49, (byte) 51, (byte) 53, (byte) 55, (byte) 57, (byte) 59, (byte) 61, (byte) 63, (byte) 65, (byte) 67, (byte) 69, (byte) 71, (byte) 73, (byte) 75, (byte) 77, (byte) 79, (byte) 81, (byte) 83, (byte) 85, (byte) 87, (byte) 89, (byte) 91, (byte) 93, (byte) 95, (byte) 97, (byte) 99, (byte) 101, (byte) 103, (byte) 105, (byte) 107, (byte) 109, (byte) 111, (byte) 113, (byte) 115, (byte) 117, (byte) 119, (byte) 121, (byte) 123, (byte) 125, (byte) 127, (byte) 129, (byte) 131, (byte) 133, (byte) 135, (byte) 137, (byte) 139, (byte) 141, (byte) 143, (byte) 145, (byte) 147, (byte) 149, (byte) 151, (byte) 153, (byte) 155, (byte) 157, (byte) 159, (byte) 161, (byte) 163, (byte) 165, (byte) 167, (byte) 169, (byte) 171, (byte) 173, (byte) 175, (byte) 177, (byte) 179, (byte) 181, (byte) 183, (byte) 185, (byte) 187, (byte) 189, (byte) 191, (byte) 193, (byte) 195, (byte) 197, (byte) 199, (byte) 201, (byte) 203, (byte) 205, (byte) 207, (byte) 209, (byte) 211, (byte) 213, (byte) 215, (byte) 217, (byte) 219, (byte) 221, (byte) 223, (byte) 225, (byte) 227, (byte) 229, (byte) 231, (byte) 233, (byte) 235, (byte) 237, (byte) 239, (byte) 241, (byte) 243, (byte) 245, (byte) 247, (byte) 249, (byte) 251, (byte) 253, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255, (byte) 255};
  128.     private final static byte[] GYR_B = new byte[] {(byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0, (byte) 0};
  129.     private final static LUT GYR_LUT = new LUT(GYR_R, GYR_G, GYR_B);
  130.    
  131.     //Variables for max Hough score search
  132.     private int maxHoughArray[][]; //Matrix to store hough scores, radii, and points from multi-threaded max search
  133.     private int ithread;
  134.     //private int totalTime = 0; //Variable to test beenfits of multithreading
  135.     // </editor-fold>
  136.    
  137.    
  138.     //Import values from GUI class before starting the analysis thread
  139.     public void setParameters(int radiusMin, int radiusMax, int radiusInc, int minCircles, int maxCircles, double thresholdRatio, int resolution, double ratio, int searchBand,
  140.                 int searchRadius, boolean reduce, boolean local, boolean houghSeries, boolean showCircles, boolean showID, boolean showScores, boolean results, boolean isGUI){
  141.        
  142.         this.radiusMin = radiusMin;
  143.         this.radiusMax = radiusMax;
  144.         this.radiusInc = radiusInc;
  145.         this.minCircles = minCircles;
  146.         this.maxCircles = maxCircles;
  147.         this.thresholdRatio = thresholdRatio;
  148.         this.resolution = resolution;
  149.         this.ratio = ratio;
  150.         this.searchBand = searchBand;
  151.         this.searchRadius = searchRadius;
  152.         this.reduce = reduce;
  153.         this.local = local;
  154.         this.houghSeries = houghSeries;
  155.         this.showCircles = showCircles;
  156.         this.showID = showID;
  157.         this.showScores = showScores;
  158.         this.results = results;
  159.         this.isGUI = isGUI;
  160.     }
  161.    
  162.     @Override
  163.     //Start the Hough transform on a separate thread from the GUI
  164.     protected Integer doInBackground() throws Exception{
  165.         startTransform();
  166.         return 100;
  167.      }
  168.    
  169.     public void startTransform(){
  170.         //Initialize the results table
  171.         rt = Analyzer.getResultsTable();
  172.         if (rt == null) {
  173.                 rt = new ResultsTable();
  174.                 Analyzer.setResultsTable(rt);
  175.         }
  176.        
  177.         //Initialize variables
  178.         nCircles = 0;
  179.         circleIDcounter = 0;
  180.        
  181.         //Calculate Hough parameters
  182.         depth = ((radiusMax-radiusMin)/radiusInc)+1;
  183.        
  184.         //If radiusInc is not a divisor of radiusMax-radiusMin, return error
  185.         if((radiusMax-radiusMin)%radiusInc != 0){
  186.             IJ.showMessage("Error: Radius increment must be a divisor of maximum radius - minimum radius.");
  187.             IJ.showMessage("radiusMin=" + radiusMin + ", radiusMax=" + radiusMax + ", radiusInc=" + radiusInc);
  188.             return;          
  189.         }
  190.        
  191.         //Build the transform LUT (all necessary translations in Cartesian coordinates)
  192.         //NOTE: This step must precede the calculatation of the threshold, as it can change the resolution from the original input value
  193.         lutSize = buildLookUpTable();
  194.          
  195.         //Calculate the threshold based off the the actual resolution
  196.         threshold = (int) Math.round(thresholdRatio * resolution);
  197.        
  198.         //If the threshold rounds to 0, set threshold to 1 as 0 is nonsensical
  199.         if(threshold <= 0) {
  200.             threshold = 1;
  201.         }        
  202.         // <editor-fold desc="Send arguments to record">
  203.         //If the macro was being recorded, return the set values to the recorder
  204.         if (Recorder.record){
  205.             String Command  = "run(\"Hough Circle Transform\",\"minRadius=" + radiusMin + ", maxRadius=" + radiusMax + ", inc=" + radiusInc +
  206.                     ", minCircles=" + minCircles + ", maxCircles=" + maxCircles + ", threshold=" + thresholdRatio + ", resolution=" + resolution +
  207.                      ", ratio=" + ratio + ", bandwidth=" + searchBand + ", local_radius=" + searchRadius + ", ";
  208.             if(reduce) Command += " reduce";
  209.             if(local) Command += " local_search";
  210.             if(houghSeries) Command += " show_raw";
  211.             if(showCircles) Command += " show_mask";
  212.             if(showID) Command += " show_centroids";
  213.             if(showScores) Command += " show_scores";
  214.             if(results) Command += " results_table";
  215.             Command += "\");\r\n";
  216.             Recorder.recordString(Command);
  217.         }
  218.         // </editor-fold>
  219.  
  220.         //Create an ImagePlus instance of the currently active image
  221.         imp = WindowManager.getCurrentImage();
  222.        
  223.         //Import the ImagePlus as a stack
  224.         ImageStack stack = imp.getStack();
  225.        
  226.         //Get pixel dimensions and units
  227.         pixelCal = imp.getCalibration();
  228.         pixelUnits = pixelCal.getUnits();
  229.         pixelDimensions = pixelCal.pixelWidth;
  230.         imageDimensions = imp.getDimensions(); //(width, height, nChannels, nSlices, nFrames)
  231.        
  232.         //If the stack has frames, get the time units and step size
  233.         if(imageDimensions[2] == 1 & imageDimensions[3] == 1 & imageDimensions[4] > 1){
  234.             timeUnits = pixelCal.getTimeUnit();
  235.             timeDimension = pixelCal.frameInterval;
  236.            
  237.             //If the time dimension is zero, then use slice # instead
  238.             if(timeDimension == 0){
  239.                 timeUnits = "slice #";
  240.                 timeDimension = 1;
  241.             }
  242.            
  243.             //Put frames into slice dimension to keep things constant
  244.             imp.setDimensions(1, imageDimensions[4], 1);
  245.         }
  246.        
  247.         //If the stack has Z-steps then use slice # as units
  248.         else if(imageDimensions[2] == 1 & imageDimensions[3] >= 1 & imageDimensions[4] == 1){
  249.             timeUnits = "slice #";
  250.             timeDimension = 1;  
  251.         }
  252.  
  253.         //If the stack is a hyper-stack, then abort analysis
  254.         else{
  255.             IJ.showMessage("Error: This stack is a hyperstack.  Please convert to stack before processing.");
  256.             IJ.showMessage("nChannels=" + imageDimensions[2] + ", nSlices=" + imageDimensions[3] + ", nFrames=" + imageDimensions[4]);
  257.             return;
  258.         }
  259.        
  260.         //Get the ROI dimensions - no ROI = full image
  261.         r = stack.getRoi();
  262.         offx = r.x;
  263.         offy = r.y;
  264.         width = r.width;
  265.         height = r.height;
  266.         fullWidth = stack.getWidth();
  267.         fullHeight = stack.getHeight();      
  268.  
  269.         //Convert the stack to float (allows 8, 16, and 32 stacks to all be processed as one type)
  270.         ImageStack stackCopy = stack.duplicate();
  271.         ImageStack floatStack = stackCopy.convertToFloat();
  272.  
  273.         if(houghSeries){
  274.             //Frames is transform dimension
  275.             houghStack = new ImageStack(width, height, depth*imp.getNSlices());
  276.         }
  277.         if(showCircles){
  278.             circleStack = new ImageStack(width, height, imp.getNSlices());
  279.         }
  280.         if(showID){
  281.             idStack = new ImageStack(width, height, imp.getNSlices());
  282.         }
  283.         if(showScores){
  284.             scoreStack = new ImageStack(width, height, imp.getNSlices());
  285.         }
  286.  
  287.         //See if input is stack
  288.         stackSlices = imp.getStackSize();
  289.         if(stackSlices > 1) isStack = true;
  290.        
  291.         //Build arrays for storing the circle result parameters
  292.         houghScores = new int [maxCircles];
  293.         centerRadii = new int [maxCircles];
  294.         centerPoint = new Point[maxCircles];
  295.         circleID = new int [maxCircles];
  296.         Arrays.fill(circleID, -1); //Flag indeces as -1 = no circle found for this index
  297.         Arrays.fill(centerRadii, -1);
  298.         Arrays.fill(centerPoint, new Point(-1,-1));
  299.         Arrays.fill(houghScores, -1);
  300.        
  301.         //Start timer to delay log outputs
  302.         long startTime = System.currentTimeMillis();
  303.        
  304.         //Retrieve the current slice in the stack as an ImageProcessor
  305.         for(int slice = 1; slice <= stackSlices; slice++){
  306.             //Store and then reset nCircles
  307.             nCirlcesPrev = nCircles;
  308.             nCircles = 0;
  309.            
  310.             //Set houghMaximum to -1 so that it can be determined whether the maximum has already been found
  311.             maxHough = -1;
  312.            
  313.             //Initialize local search to false to record whether a local search was performed
  314.             localSearch = false;
  315.          
  316.             //Show tranform status if there is more than one slice
  317.             if(isStack){
  318.                if((System.currentTimeMillis() - startTime) > GUI_UPDATE_DELAY){
  319.                    if(isGUI){
  320.                        //Update GUI progress bar
  321.                        publish("Processing frame: " + slice + " of " + stackSlices + "");
  322.                        setProgress(Math.round(100*slice/stackSlices));
  323.                    }
  324.                    
  325.                    //Update IJ status bar
  326.                    IJ.showStatus("Hough tranform - processing frame: " + slice + " of " + stackSlices + "");
  327.                    IJ.showProgress(slice, stackSlices);
  328.  
  329.                    //Reset timer
  330.                    startTime = System.currentTimeMillis();
  331.                }
  332.             }
  333.  
  334.             ImageProcessor ip = floatStack.getProcessor(slice);
  335.             imageValues = (float[]) ip.getPixels();
  336.  
  337.             // <editor-fold desc="Local search">
  338.             //If the number of circles is greater than min circles - speed up the search by looking locally for each next circle
  339.             //Multithread by giving each core a subset of the total number of circles
  340.             if(local & slice > 1){
  341.                 //If there are still a sufficient number of circles, perform an exclusively local search
  342.                 if(nCirlcesPrev>=minCircles){
  343.                     method = "Local";
  344.                    
  345.                     //Set local search to true
  346.                     localSearch = true;
  347.                     startLocalTransform();
  348.                    
  349.                     //Rebuild the local searches into a full transform if the Hough Series is desired
  350.                     if(houghSeries){
  351.                         convertLocaltoFullHough(slice, houghStack);
  352.  
  353.                     }
  354.                 }
  355.                
  356.                 //Otherwise, if there still are valid circles, perform the full Hough transform, but perform a local search for the valid circles first before continuing onto a full search
  357.                 else if(nCirlcesPrev > 0){
  358.                     method = "Partial Local";
  359.                     houghTransform();
  360.                     if (houghSeries){
  361.                         //Create the hyperstach to put into the result if needed
  362.                         HoughSpaceSeries(slice, houghStack);
  363.                     }
  364.                     startPartialLocalSearch();
  365.                     getCenterPoints();
  366.                 }
  367.                  //If there are no circles, and minCircles is 0, then return blank results
  368.                 else if(nCirlcesPrev == 0 & minCircles == 0){
  369.                     method = "N/A";
  370.                     //Create an empty Hough array
  371.                     houghValues = new int[width][height][depth];
  372.                     if (houghSeries){
  373.                         //Create the hyperstach to put into the result if needed
  374.                         HoughSpaceSeries(slice, houghStack);
  375.                     }
  376.                    
  377.                     //Set maxHough to an arbitrary value to flag that there is no need to find this value
  378.                     maxHough = 1;
  379.                 }
  380.  
  381.                 //If no circles are found, then revert to a full Hough
  382.                 else{
  383.                     method = "Full";                  
  384.                     houghTransform();
  385.                    
  386.                     if (houghSeries){
  387.                         //Create the hyperstach to put into the result if needed
  388.                         HoughSpaceSeries(slice, houghStack);
  389.                     }
  390.                     // Mark the center of the found circles in a new image(always do this as center points are necessary for local search)                  
  391.                     getCenterPoints();
  392.                 }                
  393.             }
  394.             //Otherwise, perform the full transform
  395.             else{
  396.                 method = "Full";
  397.                 houghTransform();            
  398.                 if (houghSeries){
  399.                     //Create the hyperstach to put into the result if needed
  400.                     HoughSpaceSeries(slice, houghStack);
  401.                 }
  402.                 // Mark the center of the found circles in a new image if user wants to find centers              
  403.                 if(showCircles || showID || showScores || results || local) getCenterPoints();
  404.                
  405.             }
  406.              // </editor-fold>
  407.  
  408.             // Create image View for Marked Circles.
  409.             if(showCircles) drawCircles(slice, circleStack, width, height, offx, offy, fullWidth);
  410.            
  411.             //Create map of centroids where the intensity is the radius
  412.             if (showID || showScores) drawFilledCircles(slice, idStack, scoreStack);
  413.            
  414.             //Export measurements to the results table
  415.             if (results) resultsTable(slice);
  416.            
  417.         }
  418. //startTestTime = System.currentTimeMillis();
  419. //totalTime += System.currentTimeMillis()-startTestTime;
  420. //IJ.log("" + totalTime);
  421.         //Draw the resulting stacks
  422.          if(houghSeries){
  423.             houghPlus = new ImagePlus("Hough Transform Series", houghStack);
  424.             if(depth>1) houghPlus = HyperStackConverter.toHyperStack(houghPlus, 1, depth, imp.getNSlices(), "default", "grayscale");
  425.             houghPlus.show();
  426.          }
  427.          if(showCircles){
  428.             ImagePlus Circle_Map = new ImagePlus("Centroid overlay", circleStack);
  429.             Circle_Map.show();
  430.            
  431.             //If orignal stack was movie, convert this stack to movie
  432.             if(imageDimensions[4] > 1){          
  433.                 Circle_Map.setDimensions(1, 1, imageDimensions[4]);
  434.             }            
  435.             Circle_Map.setCalibration(pixelCal);
  436.          }
  437.          if(showID){
  438.             ImagePlus ID_Map = new ImagePlus("Centroid map", idStack);
  439.             ID_Map.show();
  440.             ID_Map.setLut(RAND_LUT);
  441.             ID_Map.setDisplayRange(0,circleIDcounter);
  442.            
  443.             //If orignal stack was movie, convert this stack to movie
  444.             if(imageDimensions[4] > 1){          
  445.                 ID_Map.setDimensions(1, 1, imageDimensions[4]);
  446.             }            
  447.             ID_Map.setCalibration(pixelCal);
  448.          }
  449.          if(showScores){
  450.              ImagePlus Score_Map = new ImagePlus("Score map", scoreStack);
  451.              Score_Map.show();
  452.              Score_Map.setLut(GYR_LUT);
  453.              Score_Map.setDisplayRange(thresholdRatio, 1D);
  454.              
  455.             //If orignal stack was movie, convert this stack to movie
  456.             if(imageDimensions[4] > 1){          
  457.                 Score_Map.setDimensions(1, 1, imageDimensions[4]);
  458.             }            
  459.             Score_Map.setCalibration(pixelCal);
  460.          }
  461.          if(results) rt.show("Results");
  462.          
  463.          //Put slices back into frame dimension if necessary
  464.          if(imageDimensions[4] > 1){          
  465.             imp.setDimensions(1, 1, imageDimensions[4]);
  466.          }
  467.          
  468.          IJ.showProgress(0);
  469.     }
  470.    
  471.     //OPTMIZED - cancellable
  472.     private void startLocalTransform(){
  473.          //Initialize an array to store the local Hough transform from each thread
  474.         localHoughValues = new int [nCirlcesPrev][2*searchRadius + 1][2*searchRadius + 1][2*searchBand + 1];
  475.        
  476.         //Initialize an array to store the final dimensions and location of each local Hough Transform
  477.         localHoughParameters = new int [nCirlcesPrev][9];
  478.  
  479.         //Setup two processing pipelines to multithread only if there are more than 1 circles
  480.         if(nCirlcesPrev < 2){
  481.             localHoughTransform(0);
  482.             localGetCenterPoint(0); //Make sure this returns -1 to index if no circle was found
  483.             maxHough = houghScores[0];
  484.         }
  485.        
  486.         //Otherwise, multithread the local search
  487.         else{
  488.             //Build an array to store the result from each thread
  489.             final Thread[] threads = newThreadArray();
  490.  
  491.             //Create an atomic integer counter that each thread can use to count through the radii
  492.             final AtomicInteger ai = new AtomicInteger(0);
  493.  
  494.             //Build a thread for as many CPUs as are available to the JVM
  495.             for (ithread = 0; ithread < threads.length; ithread++) {    
  496.  
  497.                 // Concurrently run in as many threads as CPUs  
  498.                 threads[ithread] = new Thread() {  
  499.  
  500.                     { setPriority(Thread.NORM_PRIORITY); }  
  501.  
  502.                     @Override
  503.                     public void run() {  
  504.  
  505.                         //Divide the task so that each core works on a subset of circles
  506.                         for(int circleNum = ai.getAndAdd(1); circleNum < nCirlcesPrev; circleNum = ai.getAndAdd(1)){
  507.                             //Check for interrupt
  508.                             if(cancelThread) return;
  509.                            
  510.                             localHoughTransform(circleNum);
  511.                             localGetCenterPoint(circleNum); //Make sure this returns -1 to index if no circle was found
  512.                         }
  513.                     }              
  514.                 };  
  515.             }    
  516.             startAndJoin(threads);
  517.            
  518.             //Retrieve the maximum hough value if the raw Hough is desired
  519.             if(houghSeries){
  520.                 for(int i=0; i<nCirlcesPrev; i++){
  521.                     if(houghScores[i] > maxHough) maxHough = houghScores[i];
  522.                 }
  523.             }
  524.         }
  525.  
  526.         //Flag move all circles to the starting indexes, so there are no gaps between circles (i.e. if there were 3 circles, then they occupy indeces 0-2)
  527.         collapseLocalResult();
  528.     }
  529.    
  530.     //OPTIMIZED - cancellable
  531.     private void startPartialLocalSearch(){
  532.         //Build an array to store the result from each thread
  533.         final Thread[] threads = newThreadArray();
  534.  
  535.         //Create an atomic integer counter that each thread can use to count through the radii
  536.         final AtomicInteger ai = new AtomicInteger(0);  
  537.  
  538.         //Build a thread for as many CPUs as are available to the JVM
  539.         for (ithread = 0; ithread < threads.length; ithread++) {    
  540.  
  541.             // Concurrently run in as many threads as CPUs  
  542.             threads[ithread] = new Thread() {  
  543.  
  544.                 { setPriority(Thread.NORM_PRIORITY); }  
  545.  
  546.                 @Override
  547.                 public void run() {  
  548.  
  549.                     //Divide the task so that each core works on a subset of circles
  550.                     for(int circleNum = ai.getAndAdd(1); circleNum < nCirlcesPrev; circleNum = ai.getAndAdd(1)){
  551.                         getLocalCenterPoint2(circleNum);
  552.                        
  553.                         //Check for interrupt
  554.                         if(cancelThread) return;
  555.                     }
  556.                 }              
  557.             };  
  558.         }    
  559.         startAndJoin(threads);
  560.        
  561.         //Clear out the remainder of the previous circles in the circle information arrays
  562.         for(int a = nCirlcesPrev; a<circleID.length; a++){
  563.             circleID[a] = -1;
  564.             centerPoint[a] = new Point(-1,-1);
  565.             centerRadii[a] = -1;
  566.             houghScores[a] = -1;
  567.         }
  568.        
  569.         //Flag move all circles to the starting indexes, so there are no gaps between circles (i.e. if there were 3 circles, then they occupy indeces 0-2)
  570.         collapseLocalResult();
  571.     }
  572.    
  573.     //OPTMIZED - not time limiting - cancellable
  574.     //Build an array of the cartesion transformations necessary for each increment at each radius
  575.     private int buildLookUpTable() {      
  576.         //Build an array to store the X and Y coordinates for each angle increment (resolution) across each radius (depth)
  577.         lut = new int[2][resolution][depth];
  578.        
  579.         //Initialize a variable that will allow for measuring the maximium LUT size
  580.         int maxLUT = 0;
  581.        
  582.         //Step through all radii to be sampled
  583.         for(int radius = radiusMax; radius>=radiusMin; radius -= radiusInc) {
  584.             //Check for interrupt
  585.             if(cancelThread) return 0;
  586.            
  587.             //Index counter that also tracks the largest actual LUT size (may be <= resolution)
  588.             int i = 0;
  589.             for(int resStep = 0; resStep < resolution; resStep++) {
  590.                
  591.                 //Calcualte the angle and corresponding X and Y displacements for the specified radius
  592.                 double angle = (2D*Math.PI * (double)resStep) / (double)resolution;
  593.                 int indexR = (radius-radiusMin)/radiusInc;
  594.                 int rcos = (int)Math.round ((double)radius * Math.cos (angle));
  595.                 int rsin = (int)Math.round ((double)radius * Math.sin (angle));
  596.                
  597.                 //Test to make sure that the coordinate is a new coordinate
  598.                 //NOTE: A continuous circle is being discretized into pixels, it is possible that small angle steps for small radii circles will occupy the same pixel.
  599.                 //Since there is no point in making redundant calculations, these points are excluded from the LUT
  600.                 //NOTE: Using the minRadius as the transform cutoff results in a strong harmonic of the image forming near the min radius
  601.                 //threfore, using the max will push this harmonic outside the search range.
  602.                 if(radius == radiusMax && reduce){
  603.                     if(i == 0) {
  604.                         lut[0][i][indexR] = rcos;
  605.                         lut[1][i][indexR] = rsin;
  606.                         i++;
  607.                     }
  608.                     else if ((rcos != lut[0][i-1][indexR]) | (rsin != lut[1][i-1][indexR])){
  609.                         lut[0][i][indexR] = rcos;
  610.                         lut[1][i][indexR] = rsin;
  611.                         i++;
  612.                     }
  613.                 }
  614.                 else{
  615.                     lut[0][i][indexR] = rcos;
  616.                     lut[1][i][indexR] = rsin;
  617.                     i++;
  618.                 }
  619.             }
  620.            
  621.             //If this is the smallest radius, see how many transforms could be done, and set this as the new resolution
  622.             if(radius == radiusMax){
  623.                 maxLUT = i;
  624.                 resolution = maxLUT;            
  625.             }
  626.         }
  627.         return maxLUT;
  628.     }
  629.  
  630.     //OPTIMIZED - cancellable
  631.     private void houghTransform () {
  632.         //Update progress bar string with current task
  633.         if(isGUI) publish("Performing full Hough transform...");
  634.         IJ.showStatus("Performing full Hough transform...");
  635.        
  636.         //Build an array to store the result from each thread
  637.         final Thread[] threads = newThreadArray();
  638.        
  639.         //Create an atomic integer counter that each thread can use to count through the radii
  640.         final AtomicInteger ai = new AtomicInteger(radiusMin);
  641.         final AtomicInteger progress = new AtomicInteger(0);
  642.         final AtomicInteger lastProgress = new AtomicInteger(0);
  643.        
  644.         //Create an array to store the Hough values
  645.         houghValues = new int[width][height][depth];
  646.        
  647.         //Create a variable for storing the total progress possible (100 = complete)
  648.         //Nute depth is devided by nCPUs, therefore depth*nCPUs/nCPUs = depth
  649.         double totalProgress = height*depth/100;
  650.        
  651.         //Build a thread for as many CPUs as are available to the JVM
  652.         for (ithread = 0; ithread < threads.length; ithread++) {    
  653.            
  654.             // Concurrently run in as many threads as CPUs  
  655.             threads[ithread] = new Thread() {  
  656.                          
  657.                 { setPriority(Thread.NORM_PRIORITY); }  
  658.  
  659.                 @Override
  660.                 public void run() {  
  661.  
  662.                 //Divide the radius tasks across the cores available
  663.                 int currentProgress = 0;
  664.                 for (int radius = ai.getAndAdd(radiusInc); radius <= radiusMax; radius = ai.getAndAdd(radiusInc)) {  
  665.                     int indexR=(radius-radiusMin)/radiusInc;
  666.                     //For a given radius, transform each pixel in a circle, and add-up the votes
  667.                     for(int y = 1; y < height-1; y++) {
  668.                         //Increment the progress counter, and submit the current progress status
  669.                         progress.getAndAdd(1);
  670.                        
  671.                         //Calculate the current progress value
  672.                         currentProgress = Math.round((float) (progress.get()/totalProgress));
  673.  
  674.                         //There is a significant time penalty for progress updates, so only update if needed
  675.                         if(currentProgress > lastProgress.get()){ //7.8s with if, 8.7s without if, 7.8s with no progress update, 8.7s with delay between GUI updates
  676.                             if(isGUI && currentProgress <= 100) setProgress(currentProgress);
  677.                             IJ.showProgress(currentProgress, 100);
  678.                             lastProgress.set(currentProgress);
  679.                         }
  680.                        
  681.                         //Check for interrupt
  682.                         if(cancelThread) return;
  683.                        
  684.                         for(int x = 1; x < width-1; x++) {
  685.                                 if( imageValues[(x+offx)+(y+offy)*fullWidth] != 0 )  {// Edge pixel found                                    
  686.                                     for(int i = 0; i < lutSize; i++) {
  687.                                         int a = x + lut[1][i][indexR];
  688.                                         int b = y + lut[0][i][indexR];
  689.                                         if((b >= 0) & (b < height) & (a >= 0) & (a < width)) {
  690.                                             houghValues[a][b][indexR] += 1;
  691.                                         }
  692.                                     }
  693.                                 }
  694.                             }
  695.                         }
  696.                     }
  697.                 }  
  698.             };  
  699.         }    
  700.         startAndJoin(threads);      
  701.     }
  702.    
  703.     //The local Hough is an inversion of the inertial reference frame used in the full Hough
  704.     //In the full Hough the image is seach for pixels with a value > 1, and if this is true
  705.     //then the pixel is projected in a circle about that point.
  706.    
  707.     //OPTMIZED - cancellable
  708.     //To reduce the necessary transform space,
  709.     private void localHoughTransform (int index) {
  710.         //Initialize local search variables
  711.         int startWidth = centerPoint[index].x - searchRadius;
  712.         if (startWidth < 1) startWidth = 1; //Keep search area inside the image
  713.         int endWidth = centerPoint[index].x + searchRadius;
  714.         if (endWidth > width) endWidth = width; //Keep search area inside the image
  715.         int lwidth = endWidth - startWidth + 1;
  716.        
  717.         int startHeight = centerPoint[index].y - searchRadius;
  718.         if (startHeight < 1) startHeight = 1; //Keep search area inside the image
  719.         int endHeight = centerPoint[index].y + searchRadius;
  720.         if (endHeight > height) endHeight = height; //Keep search area inside the image
  721.         int lheight = endHeight - startHeight + 1;
  722.        
  723.         int lradiusMin = centerRadii[index] - searchBand;
  724.         if(lradiusMin < radiusMin) lradiusMin = radiusMin;
  725.         int lradiusMax = centerRadii[index] + searchBand;
  726.         if(lradiusMax > radiusMax) lradiusMax = radiusMax;
  727.         int ldepth = ((lradiusMax-lradiusMin)/radiusInc)+1;
  728.        
  729.         //Store local Hough parameters into array
  730.         localHoughParameters[index][0] = startWidth;
  731.         localHoughParameters[index][1] = endWidth;
  732.         localHoughParameters[index][2] = lwidth;
  733.         localHoughParameters[index][3] = startHeight;
  734.         localHoughParameters[index][4] = endHeight;
  735.         localHoughParameters[index][5] = lheight;
  736.         localHoughParameters[index][6] = lradiusMin;
  737.         localHoughParameters[index][7] = lradiusMax;
  738.         localHoughParameters[index][8] = ldepth;
  739.        
  740.         //Divide the radius tasks across the cores available
  741.         for (int radius = lradiusMin; radius <= lradiusMax; radius += radiusInc) {
  742.             int indexR=(radius-lradiusMin)/radiusInc;
  743.             //For a given radius, transform each pixel in a circle, and add-up the votes
  744.             for(int y = startHeight; y < endHeight-1; y++) {
  745.                
  746.                 //Check for interrupt
  747.                 if(cancelThread) return;
  748.                
  749.                 for(int x = startWidth; x < endWidth-1; x++) {
  750.                     for(int i = 0; i < lutSize; i++){
  751.                         int a = x + offx + lut[1][i][indexR+((lradiusMin-radiusMin)/radiusInc)];
  752.                         int b = y + offy + lut[0][i][indexR+((lradiusMin-radiusMin)/radiusInc)];
  753.                        
  754.                         //Check to make sure pixel is within the image
  755.                         if((a>1) & (a<fullWidth-1) & (b>1) & (b<fullHeight-1)){
  756.                             //See if the pixel has an intensity > 1, if so, add 1 to the vote
  757.                             if(imageValues[a + (b*fullWidth)] > 0 ){
  758.                                 localHoughValues[index][x-startWidth][y-startHeight][indexR] += 1;
  759.                             }
  760.                         }
  761.                     }
  762.                 }
  763.             }
  764.         }
  765.     }  
  766.  
  767.     //OPTMIZED - cancellable
  768.     //Find the largest Hough pixel in the 3D Hough transform array to scale the 8-bit conversion
  769.     private void houghMaximum () {
  770.         //long startTime = System.currentTimeMillis(); //1337ms without multi, 319ms with multi, 175ms by writing to private variable per thread
  771.         //Build an array to store the result from each thread
  772.         final Thread[] threads = newThreadArray();
  773.        
  774.         //Build an array to store the max from each thread
  775.         maxHoughArray = new int[threads.length][4];
  776.        
  777.         //Create an atomic integer counter that each thread can use to count through the radii
  778.         final AtomicInteger ai = new AtomicInteger(0);
  779.         final AtomicInteger progress = new AtomicInteger(0);
  780.         final AtomicInteger lastProgress = new AtomicInteger(0);
  781.                  
  782.         //Create an integer for indexing the results array (one result per thread
  783.         final AtomicInteger Index = new AtomicInteger(0);
  784.        
  785.         //Create a variable for storing the total progress possible (100 = complete)
  786.         //Nute depth is devided by nCPUs, therefore depth*nCPUs/nCPUs = depth
  787.         double totalProgress = height*depth/100;
  788.  
  789.         maxHough = 0;
  790.        
  791.         //Build a thread for as many CPUs as are available to the JVM
  792.         for (ithread = 0; ithread < threads.length; ithread++) {    
  793.            
  794.             // Concurrently run in as many threads as CPUs  
  795.             threads[ithread] = new Thread() {  
  796.                          
  797.                 { setPriority(Thread.NORM_PRIORITY); }  
  798.  
  799.                 //Search for the largest score with each thread
  800.                 @Override
  801.                 public void run() {
  802.                     int maxHoughThread = -1;
  803.                     int maxRadiusThread = -1;
  804.                     int currentProgress = 0;
  805.                     Point maxPointThread = new Point (-1,-1);
  806.                     for(int a=ai.getAndIncrement(); a<depth; a=ai.getAndIncrement()){
  807.                         for(int j = 0; j < height; j++) {
  808.                             //Increment the progress counter, and submit the current progress status
  809.                             progress.getAndAdd(1);
  810.                            
  811.                             //Gui updates can be time intensive, so only update at fixed time intervals
  812.                             currentProgress = Math.round((float) (progress.get()/totalProgress));
  813.  
  814.                             //There is a significant time penalty for progress updates, so only update if needed
  815.                             if(currentProgress > lastProgress.get() & currentProgress >= 0 & currentProgress <= 100){ //7.8s with if, 8.7s without if, 7.8s with no progress update, 8.7s with delay between GUI updates
  816.                                 if(isGUI) setProgress(currentProgress);
  817.                                 IJ.showProgress(currentProgress, 100);
  818.                                 lastProgress.set(currentProgress);
  819.                             }
  820.                            
  821.                             //Check for interrupt
  822.                             if(cancelThread) return;
  823.                            
  824.                             for(int k = 0; k < width; k++){
  825.                                 if(houghValues[k][j][a] > maxHoughThread) {
  826.                                     maxHoughThread = houghValues[k][j][a];
  827.                                     maxPointThread = new Point(k,j);
  828.                                     maxRadiusThread = a*radiusInc + radiusMin;
  829.                                 }
  830.                             }
  831.                         }
  832.                     }
  833.                     //Have each thread report the score to a common array
  834.                     maxHoughArray[Index.getAndIncrement()] = new int[]{maxHoughThread, maxRadiusThread, maxPointThread.x, maxPointThread.y};                    
  835.                 }
  836.             };
  837.         }
  838.         startAndJoin(threads);
  839.        
  840.         //Search common array for highest score
  841.         for (int[] maxHoughArray1 : maxHoughArray) {
  842.             if (maxHoughArray1[0] > maxHough) {
  843.                 maxHough = maxHoughArray1[0];
  844.                 maxRadius = maxHoughArray1[1];
  845.                 maxPoint = new Point((int) maxHoughArray1[2], (int) maxHoughArray1[3]);    
  846.             }
  847.         }
  848.     }
  849.    
  850.     //OPTMIZED - cancellable
  851.     //Create a Hough stack series using the local transforms
  852.     private void convertLocaltoFullHough(int slice, ImageStack houghStack){
  853.        
  854.        
  855.         int startFrame = ((slice-1)*depth);
  856.  
  857.         //Create an array to store the Hough values
  858.         byte localHoughPixels[][] = new byte[depth][width*height];
  859.        
  860.         //NOTE: Only do multithreading if the problem is sufficiently complex
  861.         if(localHoughParameters.length < 2 & localHoughParameters.length*searchBand*searchRadius*searchRadius < 15000){
  862.             //Extract the local search parameters
  863.             //Divide the task so that each core works on a subset of circles
  864.             for(int circleNum = 0; circleNum < localHoughParameters.length; circleNum++){
  865.                 //Extract the local search parameters
  866.                 int startWidth = localHoughParameters[circleNum][0];
  867.                 int endWidth = localHoughParameters[circleNum][1];
  868.                 int lwidth = localHoughParameters[circleNum][2];
  869.                 int startHeight = localHoughParameters[circleNum][3];
  870.                 int endHeight = localHoughParameters[circleNum][4];
  871.                 int lheight = localHoughParameters[circleNum][5];
  872.                 int lradiusMin = localHoughParameters[circleNum][6];
  873.                 int lradiusMax = localHoughParameters[circleNum][7];
  874.                 int ldepth = localHoughParameters[circleNum][8];
  875.                 for(int h=0; h<lheight; h++){
  876.                     //Check for interrupt
  877.                     if(cancelThread) return;
  878.                    
  879.                     for(int w=0; w<lwidth; w++){
  880.                         for(int i=0; i<ldepth; i++){
  881.                            localHoughPixels[i+(lradiusMin-radiusMin)/radiusInc][(w+startWidth)+(h+startHeight)*width] = (byte) Math.round ((localHoughValues[circleNum][w][h][i] * 255D) / maxHough);
  882.                         }
  883.                     }
  884.                 }
  885.             }
  886.         }
  887.        
  888.         //Otherwise, multithread the depositing of pixels
  889.         else{
  890.             //Build an array to store the result from each thread
  891.             final Thread[] threads = newThreadArray();
  892.  
  893.             //Create an atomic integer counter that each thread can use to count through the radii
  894.             final AtomicInteger ai = new AtomicInteger(0);  
  895.  
  896.             //Build a thread for as many CPUs as are available to the JVM
  897.             for (ithread = 0; ithread < threads.length; ithread++) {    
  898.  
  899.                 // Concurrently run in as many threads as CPUs  
  900.                 threads[ithread] = new Thread() {  
  901.  
  902.                     { setPriority(Thread.NORM_PRIORITY); }  
  903.  
  904.                     @Override
  905.                     public void run() {  
  906.  
  907.                         //Divide the task so that each core works on a subset of circles
  908.                         for(int circleNum = ai.getAndAdd(1); circleNum < localHoughParameters.length; circleNum = ai.getAndAdd(1)){
  909.                             //Extract the local search parameters
  910.                             int startWidth = localHoughParameters[circleNum][0];
  911.                             int endWidth = localHoughParameters[circleNum][1];
  912.                             int lwidth = localHoughParameters[circleNum][2];
  913.                             int startHeight = localHoughParameters[circleNum][3];
  914.                             int endHeight = localHoughParameters[circleNum][4];
  915.                             int lheight = localHoughParameters[circleNum][5];
  916.                             int lradiusMin = localHoughParameters[circleNum][6];
  917.                             int lradiusMax = localHoughParameters[circleNum][7];
  918.                             int ldepth = localHoughParameters[circleNum][8];
  919.                             int index = 0;
  920.                             int maxIndex = width * height;
  921.                             for(int h=0; h<lheight; h++){
  922.                                 //Check for interrupt
  923.                                 if(cancelThread) return;
  924.                                
  925.                                 for(int w=0; w<lwidth; w++){
  926.                                     for(int i=0; i<ldepth; i++){
  927.                                        index = (w+startWidth)+(h+startHeight)*width;
  928.                                        if(index < maxIndex) localHoughPixels[i+(lradiusMin-radiusMin)/radiusInc][index] = (byte) Math.round ((localHoughValues[circleNum][w][h][i] * 255D) / maxHough);
  929.                                     }
  930.                                 }
  931.                             }
  932.                         }
  933.                     }              
  934.                 };  
  935.             }    
  936.             startAndJoin(threads);
  937.         }
  938.        
  939.         //Deposit the array into the Hough Series stack
  940.         //Not time limiting, even at fairly high resolutions
  941.         for(int radius = radiusMin; radius<=radiusMax; radius += radiusInc) {
  942.             //Calculate the corresponding index
  943.             int houghIndex = (radius-radiusMin)/radiusInc;
  944.            
  945.             //Deposit the array image into the corresponding slice in the stack
  946.             houghStack.setPixels(localHoughPixels[houghIndex], houghIndex+1+startFrame);
  947.  
  948.             //Give the current slice the appropriate radius label
  949.             houghStack.setSliceLabel("Hough Space [r="+radius+", resolution="+resolution+"]", houghIndex+1+startFrame);
  950.         }
  951.        
  952.     }
  953.    
  954.     //OPTIMIZED - cancellable
  955.     //Add transform series data to the hyperstack
  956.     private void HoughSpaceSeries(int slice, ImageStack houghStack){  
  957.         //If the maximum Hough value has not yet been assigned, search the whole Hough transform for the maximum value
  958.         if(maxHough == -1){
  959.             houghMaximum();  
  960.         }
  961.  
  962.         //Create an array to store the Hough values
  963.         byte localHoughPixels[][] = new byte[depth][width*height];        
  964.         int startFrame = ((slice-1)*depth);
  965.  
  966.         //If a full transform was done, calcualte the Hough Pixels
  967.         if(minCircles > 0 | nCirlcesPrev > 0){
  968.             //Build an array to store the result from each thread
  969.             final Thread[] threads = newThreadArray();
  970.  
  971.             //Create an atomic integer counter that each thread can use to count through the radii
  972.             final AtomicInteger ai = new AtomicInteger(radiusMin);  
  973.  
  974.             //Build a thread for as many CPUs as are available to the JVM
  975.             for (ithread = 0; ithread < threads.length; ithread++) {    
  976.  
  977.                 // Concurrently run in as many threads as CPUs  
  978.                 threads[ithread] = new Thread() {  
  979.  
  980.                     { setPriority(Thread.NORM_PRIORITY); }  
  981.  
  982.                     @Override
  983.                     public void run() {  
  984.  
  985.                     //Divide the radius tasks across the cores available
  986.                     for (int radius = ai.getAndAdd(radiusInc); radius <= radiusMax; radius = ai.getAndAdd(radiusInc)) {
  987.                         //Check for interrupt
  988.                         if(cancelThread) return;
  989.                        
  990.                         //Calculate the corresponding index
  991.                         int houghIndex = (radius-radiusMin)/radiusInc;
  992.  
  993.                         //If full tansform was performed, retrieve the pixel array for the current Hough radius image
  994.                         createHoughPixels(localHoughPixels[houghIndex], houghIndex);
  995.  
  996.                         //Deposit the array image into the corresponding slice in the stack
  997.                         houghStack.setPixels(localHoughPixels[houghIndex], houghIndex+1+startFrame);
  998.  
  999.                         //Give the current slice the appropriate radius label
  1000.                         houghStack.setSliceLabel("Hough Space [r="+radius+", resolution="+resolution+"]", houghIndex+1+startFrame);
  1001.                     }
  1002.                 }  
  1003.             };
  1004.             }
  1005.             startAndJoin(threads);
  1006.         }
  1007.    
  1008.         //Deposit the array into the Hough Series stack
  1009.         //Not time limiting, even at fairly high resolutions
  1010.         for(int radius = radiusMin; radius<=radiusMax; radius += radiusInc) {
  1011.             //Check for interrupt
  1012.             if(cancelThread) return;
  1013.  
  1014.             //Calculate the corresponding index
  1015.             int houghIndex = (radius-radiusMin)/radiusInc;
  1016.            
  1017.             //Deposit the array image into the corresponding slice in the stack
  1018.             houghStack.setPixels(localHoughPixels[houghIndex], houghIndex+1+startFrame);
  1019.  
  1020.             //Give the current slice the appropriate radius label
  1021.             houghStack.setSliceLabel("Hough Space [r="+radius+", resolution="+resolution+"]", houghIndex+1+startFrame);
  1022.         }
  1023.     }
  1024.    
  1025.     //OPTMIZED - cancellable
  1026.     // Convert Values in Hough Space to an 8-Bit Image Space.
  1027.     private void createHoughPixels (byte houghPixels[], int index) {
  1028.         //Rescale all the Hough values to 8-bit to create the Hough image - 47ms to complete - single threading okay
  1029.         for(int l = 0; l < height; l++) {
  1030.             //Check for interrupt
  1031.             if(cancelThread) return;
  1032.            
  1033.             for(int i = 0; i < width; i++) {
  1034.                 houghPixels[i + l * width] = (byte) Math.round ((houghValues[i][l][index] * 255D) / maxHough);
  1035.             }
  1036.  
  1037.         }
  1038.     }
  1039.  
  1040.     // Draw the circles found in the original image. - cancellable
  1041.     private void drawCircles(int slice, ImageStack circleStack, int widthROI, int heightROI, int fullWidthX, int fullWidthY, int fullWidthROI) {
  1042.                
  1043.             // Copy original input pixels into output
  1044.             // circle location display image and
  1045.             // combine with saturation at 100
  1046.             byte[] circlespixels = new byte [widthROI*heightROI];
  1047.  
  1048.             int roiaddr=0;
  1049.             for( int y = fullWidthY; y < fullWidthY+heightROI; y++) {
  1050.                     //Check for interrupt
  1051.                     if(cancelThread) return;
  1052.                     for(int x = fullWidthX; x < fullWidthX+widthROI; x++) {
  1053.                             // Copy;
  1054.                             circlespixels[roiaddr] = (byte) imageValues[x+fullWidthROI*y];
  1055.                             // Saturate
  1056.                             if(circlespixels[roiaddr] != 0 )
  1057.                                     circlespixels[roiaddr] = 100;
  1058.                             else
  1059.                                     circlespixels[roiaddr] = 0;
  1060.                             roiaddr++;
  1061.                     }
  1062.             }
  1063.             // Copy original image to the circlespixels image.
  1064.             // Changing pixels values to 100, so that the marked
  1065.             // circles appears more clear. Must be improved in
  1066.             // the future to show the resuls in a colored image.
  1067.             //for(int i = 0; i < width*height ;++i ) {
  1068.             //if(imageValues[i] != 0 )
  1069.             //if(circlespixels[i] != 0 )
  1070.             //circlespixels[i] = 100;
  1071.             //else
  1072.             //circlespixels[i] = 0;
  1073.             //}
  1074.             if(centerPoint == null) getCenterPoints();
  1075.  
  1076.             byte cor = (byte) 255;
  1077.             // Redefine these so refer to ROI coordinates exclusively
  1078.             fullWidthROI= widthROI;
  1079.             fullWidthX=0;
  1080.             fullWidthY=0;
  1081.  
  1082.             for(int l = 0; l < nCircles; l++) {
  1083.                     int i = centerPoint[l].x;
  1084.                     int j = centerPoint[l].y;
  1085.                    
  1086.                     //Check for interrupt
  1087.                     if(cancelThread) return;
  1088.                     // Draw a gray cross marking the center of each circle.
  1089.                     for( int k = -10 ; k <= 10 ; ++k ) {
  1090.                             int p = (j+k+fullWidthY)*fullWidthROI+ (i+fullWidthX);
  1091.                             if(!outOfBounds(j+k+fullWidthY,i+fullWidthX))
  1092.                                     circlespixels[(j+k+fullWidthY)*fullWidthROI+ (i+fullWidthX)] = cor;
  1093.                             if(!outOfBounds(j+fullWidthY,i+k+fullWidthX))
  1094.                                     circlespixels[(j+fullWidthY)*fullWidthROI  + (i+k+fullWidthX)] = cor;
  1095.                     }
  1096.                     for( int k = -2 ; k <= 2 ; ++k ) {
  1097.                             if(!outOfBounds(j-2+fullWidthY,i+k+fullWidthX))
  1098.                                     circlespixels[(j-2+fullWidthY)*fullWidthROI+ (i+k+fullWidthX)] = cor;
  1099.                             if(!outOfBounds(j+2+fullWidthY,i+k+fullWidthX))
  1100.                                     circlespixels[(j+2+fullWidthY)*fullWidthROI+ (i+k+fullWidthX)] = cor;
  1101.                             if(!outOfBounds(j+k+fullWidthY,i-2+fullWidthX))
  1102.                                     circlespixels[(j+k+fullWidthY)*fullWidthROI+ (i-2+fullWidthX)] = cor;
  1103.                             if(!outOfBounds(j+k+fullWidthY,i+2+fullWidthX))
  1104.                                     circlespixels[(j+k+fullWidthY)*fullWidthROI+ (i+2+fullWidthX)] = cor;
  1105.                     }
  1106.                    
  1107.                     //Draw Bresenham circle
  1108.                     int x=centerRadii[l];
  1109.                     int y=0;
  1110.                     int err = 0;
  1111.                     while (x >= y) {
  1112.                         if(!outOfBounds(j+y+fullWidthY,i+x+fullWidthX))
  1113.                             circlespixels[(j+y+fullWidthY)*fullWidthROI+ (i+x+fullWidthX)] = cor;
  1114.                         if(!outOfBounds(j+y+fullWidthY,i-x+fullWidthX))
  1115.                             circlespixels[(j+y+fullWidthY)*fullWidthROI+ (i-x+fullWidthX)] = cor;
  1116.                         if(!outOfBounds(j-y+fullWidthY,i+x+fullWidthX))
  1117.                             circlespixels[(j-y+fullWidthY)*fullWidthROI+ (i+x+fullWidthX)] = cor;
  1118.                         if(!outOfBounds(j-y+fullWidthY,i-x+fullWidthX))
  1119.                             circlespixels[(j-y+fullWidthY)*fullWidthROI+ (i-x+fullWidthX)] = cor;
  1120.  
  1121.                         if(!outOfBounds(j+x+fullWidthY,i+y+fullWidthX))
  1122.                             circlespixels[(j+x+fullWidthY)*fullWidthROI+ (i+y+fullWidthX)] = cor;
  1123.                         if(!outOfBounds(j+x+fullWidthY,i-y+fullWidthX))
  1124.                             circlespixels[(j+x+fullWidthY)*fullWidthROI+ (i-y+fullWidthX)] = cor;                    
  1125.                         if(!outOfBounds(j-x+fullWidthY,i+y+fullWidthX))
  1126.                             circlespixels[(j-x+fullWidthY)*fullWidthROI+ (i+y+fullWidthX)] = cor;                  
  1127.                         if(!outOfBounds(j-x+fullWidthY,i-y+fullWidthX))
  1128.                             circlespixels[(j-x+fullWidthY)*fullWidthROI+ (i-y+fullWidthX)] = cor;  
  1129.                        
  1130.                         if(err <= 0){
  1131.                             y += 1;
  1132.                             err += 2*y + 1;
  1133.                         }
  1134.                         else{
  1135.                             x -= 1;
  1136.                             err -= 2*x + 1;
  1137.                         }
  1138.                     }
  1139.                    
  1140.             }
  1141.             circleStack.setPixels(circlespixels, slice);
  1142.             circleStack.setSliceLabel(nCircles + " circles found", slice);
  1143.     }
  1144.    
  1145.    
  1146.     // Draw the centroids found in the original image where intensity = radius.
  1147.     private void drawFilledCircles(int slice, ImageStack idStack, ImageStack scoreStack) {
  1148.            
  1149.         //Create arrays the same size as the images
  1150.         short[] IDpixels = new short[width*height];
  1151.         float[] scorepixels = new float[width*height];
  1152.        
  1153.         for(int l = 0; l < nCircles; l++) {
  1154.             int i = centerPoint[l].x;
  1155.             int j = centerPoint[l].y;
  1156.             int radius = centerRadii[l];
  1157.             short ID = (short) circleID[l];
  1158.             float score = (float) houghScores[l]/(float) resolution;
  1159.             int rSquared = radius*radius;
  1160.            
  1161.             for(int y=-1*radius; y<=radius; y++){
  1162.                 for(int x=-1*radius; x<=radius; x++){
  1163.                     if(x*x+y*y <= rSquared){
  1164.                         if(showID){
  1165.                             if(!outOfBounds(j+y,i+x)) IDpixels[(j+y)*width + i + x] = ID;
  1166.                         }
  1167.                         if(showScores){
  1168.                             if(!outOfBounds(j+y,i+x)) scorepixels[(j+y)*width + i + x] = score;
  1169.                         }
  1170.                     }
  1171.                 }
  1172.             }
  1173.         }
  1174.        
  1175.         if(showID){
  1176.             idStack.setPixels(IDpixels, slice);
  1177.             idStack.setSliceLabel(nCircles + " circles found", slice);
  1178.         }
  1179.         if(showScores){
  1180.             scoreStack.setPixels(scorepixels, slice);
  1181.             scoreStack.setSliceLabel(nCircles + " circles found", slice);
  1182.         }
  1183.     }
  1184.    
  1185.     //Export the results to the results table
  1186.     public void resultsTable(int frame){
  1187.         for(int a = 0; a < nCircles; a++){
  1188.             rt.incrementCounter();
  1189.             rt.addValue("ID", circleID[a]);
  1190.             rt.addValue("X (" + pixelUnits + ")" , centerPoint[a].x*pixelDimensions);
  1191.             rt.addValue("Y (" + pixelUnits + ")", centerPoint[a].y*pixelDimensions);
  1192.             rt.addValue("Radius (" + pixelUnits + ")", centerRadii[a]*pixelDimensions);
  1193.             rt.addValue("Score", ((float) houghScores[a]/resolution));
  1194.             rt.addValue("nCircles", nCircles);
  1195.             rt.addValue("Resolution", lutSize);
  1196.             rt.addValue("Frame (" + timeUnits + ")", frame*timeDimension);
  1197.             rt.addValue("Method", method);
  1198.         }
  1199.     }
  1200.    
  1201.  
  1202.     private boolean outOfBounds(int y,int x) {
  1203.         if(x >= width)
  1204.             return(true);
  1205.         if(x <= 0)
  1206.             return(true);
  1207.         if(y >= height)
  1208.             return(true);
  1209.         if(y <= 0)
  1210.             return(true);
  1211.         return(false);
  1212.     }
  1213.  
  1214.     //OPTMIZED - cancellable
  1215.     /** Search for a fixed number of circles.
  1216.     @param maxCircles The number of circles that should be found.  
  1217.     */
  1218.     private void getCenterPoints () {
  1219.    
  1220.         int countCircles = nCircles;
  1221.         maxHough = threshold;
  1222.  
  1223.         while(countCircles < maxCircles && maxHough >= threshold) {
  1224.             //Check for interrupt
  1225.             if(cancelThread) return;
  1226.            
  1227.             //Update bar string with current circle that is being searched for
  1228.             if(isGUI) publish("Searching for circles. " + countCircles + " circles found.");
  1229.             IJ.showStatus("Searching for circles. " + countCircles + " circles found.");
  1230.            
  1231.             //Search for the highest remaining Hough score in the matrix
  1232.             houghMaximum();
  1233.  
  1234.             if(maxHough>=threshold){
  1235.                 circleIDcounter++;
  1236.                 circleID[countCircles] = circleIDcounter;
  1237.                 centerPoint[countCircles] = maxPoint;
  1238.                 centerRadii[countCircles] = maxRadius;
  1239.                 houghScores[countCircles] = maxHough;
  1240.                 countCircles++;
  1241.                 clearNeighbours(maxPoint.x,maxPoint.y, maxRadius);
  1242.             }
  1243.         }
  1244.  
  1245.         nCircles = countCircles;
  1246.        
  1247.         //Clear the remainder of result arrays, since these results are from the previous round
  1248.         for(int i = nCircles; i<maxCircles; i++){
  1249.             circleID[i] = -1;
  1250.             centerPoint[i] = new Point(-1,-1);
  1251.             centerRadii[i] = -1;
  1252.             houghScores[i] = -1;
  1253.         }
  1254.     }
  1255.    
  1256.     //OPTIMIZED - cancellable
  1257.     private void getLocalCenterPoint2(int index){
  1258.         //Initialize local search variables
  1259.         //Initialize local search variables
  1260.         int startWidth = centerPoint[index].x - searchRadius;
  1261.         if (startWidth < 1) startWidth = 1; //Keep search area inside the image
  1262.         int endWidth = centerPoint[index].x + searchRadius;
  1263.         if (endWidth > width) endWidth = width; //Keep search area inside the image
  1264.        
  1265.         int startHeight = centerPoint[index].y - searchRadius;
  1266.         if (startHeight < 1) startHeight = 1; //Keep search area inside the image
  1267.         int endHeight = centerPoint[index].y + searchRadius;
  1268.         if (endHeight > height) endHeight = height; //Keep search area inside the image
  1269.        
  1270.         int lradiusMin = centerRadii[index] - searchBand;
  1271.         if(lradiusMin < radiusMin) lradiusMin = radiusMin;
  1272.         int lradiusMax = centerRadii[index] + searchBand;
  1273.         if(lradiusMax > radiusMax) lradiusMax = radiusMax;
  1274.  
  1275.        
  1276.         //Search locally for the highest hough score in the full Hough Space
  1277.         int maxLocalHough = -1;
  1278.         int maxLocalRadius = -1;
  1279.         Point maxLocalPoint = new Point (-1,-1);
  1280.         for(int radius = lradiusMin; radius<=lradiusMax; radius += radiusInc){
  1281.             int indexR=(radius-radiusMin)/radiusInc;
  1282.             for(int j = startHeight; j < endHeight; j++) {
  1283.                 //Check for interrupt
  1284.                 if(cancelThread) return;
  1285.                
  1286.                 for(int k = startWidth; k < endWidth; k++){
  1287.                     if(houghValues[k][j][indexR] > maxLocalHough) {
  1288.                         maxLocalHough = houghValues[k][j][indexR];
  1289.                         maxLocalPoint = new Point(k,j);
  1290.                         maxLocalRadius = radius;
  1291.                     }
  1292.                 }
  1293.             }
  1294.         }      
  1295.         //If the highest score is above the threshold, record the new circle
  1296.         if(maxLocalHough>=threshold){
  1297.             centerPoint[index] = maxLocalPoint;
  1298.             centerRadii[index] = maxLocalRadius;
  1299.             houghScores[index] = maxLocalHough;
  1300.             clearNeighbours(maxLocalPoint.x,maxLocalPoint.y, maxLocalRadius);
  1301.         }
  1302.         //Otherwise, record that the circle was lost
  1303.         else{
  1304.             circleID[index] = -1;
  1305.             centerPoint[index] = new Point(-1,-1);
  1306.             centerRadii[index] = -1;
  1307.             houghScores[index] = -1;
  1308.         }  
  1309.     }
  1310.    
  1311.     //OPTMIZED - cancellable
  1312.     private void localGetCenterPoint(int index){
  1313.        
  1314.         //Extract the local search parameters
  1315.         int startWidth = localHoughParameters[index][0];
  1316.         int endWidth = localHoughParameters[index][1];
  1317.         int lwidth = localHoughParameters[index][2];
  1318.         int startHeight = localHoughParameters[index][3];
  1319.         int endHeight = localHoughParameters[index][4];
  1320.         int lheight = localHoughParameters[index][5];
  1321.         int lradiusMin = localHoughParameters[index][6];
  1322.         int lradiusMax = localHoughParameters[index][7];
  1323.         int ldepth = localHoughParameters[index][8];
  1324.        
  1325.         //Search for the highest hough score in the local Hough Space
  1326.         int maxLocalHough = -1;
  1327.         int maxLocalRadius = -1;
  1328.         Point maxLocalPoint = new Point (-1,-1);
  1329.         for(int a=0; a<ldepth; a++){
  1330.             for(int j = 0; j < lheight; j++) {
  1331.                 //Check for interrupt
  1332.                 if(cancelThread) return;
  1333.                
  1334.                 for(int k = 0; k < lwidth; k++){
  1335.                     if(localHoughValues[index][k][j][a] > maxLocalHough) {
  1336.                         maxLocalHough = localHoughValues[index][k][j][a];
  1337.                         maxLocalPoint = new Point(k + startWidth,j + startHeight);
  1338.                         maxLocalRadius = a*radiusInc + lradiusMin;
  1339.                     }
  1340.                 }
  1341.             }
  1342.         }
  1343.        
  1344.         //If the highest score is above the threshold, record the new circle
  1345.         if(maxLocalHough>=threshold){            
  1346.             centerPoint[index] = maxLocalPoint;
  1347.             centerRadii[index] = maxLocalRadius;
  1348.             houghScores[index] = maxLocalHough;
  1349.         }
  1350.         //Otherwise, record that the circle was lost
  1351.         else{
  1352.             circleID[index] = -1;
  1353.             centerPoint[index] = new Point(-1,-1);
  1354.             centerRadii[index] = -1;
  1355.             houghScores[index] = -1;
  1356.         }    
  1357.     }
  1358.  
  1359.     //OPTMIZED
  1360.     private void collapseLocalResult(){
  1361.         //search for indeces containing circles (i.e. index != -1)
  1362.         int[] idIndeces = new int [circleID.length];
  1363.         int indexCounter = 0;
  1364.         for(int i=0; i<circleID.length; i++){
  1365.             //If a valid ID is found, record the index
  1366.             if(circleID[i] != -1){
  1367.                 idIndeces[indexCounter] = i;
  1368.                 indexCounter++;
  1369.             }
  1370.         }
  1371.        
  1372.         //Move all the found results to the starting indeces of the arrays if needed
  1373.         if(indexCounter < maxCircles){
  1374.             for(int i=0; i<indexCounter; i++){
  1375.                 //If index is empty, then move the result
  1376.                 if(circleID[i] == -1){
  1377.                     //Check to see index is empty
  1378.                     //Move results
  1379.                     circleID[i] = circleID[idIndeces[i]];
  1380.                     houghScores[i] = houghScores[idIndeces[i]];
  1381.                     centerRadii[i] = centerRadii[idIndeces[i]];
  1382.                     centerPoint[i] = centerPoint[idIndeces[i]];
  1383.  
  1384.                     //Clear original index
  1385.                     circleID[idIndeces[i]] = -1;
  1386.                     houghScores[idIndeces[i]] = -1;
  1387.                     centerRadii[idIndeces[i]] = -1;
  1388.                     centerPoint[idIndeces[i]] = new Point(-1,-1);
  1389.                 }
  1390.             }
  1391.         }
  1392.        
  1393.         //Update nCircles to reflect the new found number of circles
  1394.         nCircles = indexCounter;
  1395.     };
  1396.    
  1397.     //OPTIMIZED - Not time limiting, even with large circles - cancellable
  1398.     /** Clear, from the Hough Space, all the counter that are near (radius/2) a previously found circle C.    
  1399.     @param x The x coordinate of the circle C found.
  1400.     @param x The y coordinate of the circle C found.
  1401.     @param x The radius of the circle C found.
  1402.     */
  1403.     private void clearNeighbours(int x,int y, int radius) {
  1404.         // The following code just clean the points around the center of the circle found.
  1405.         int radiusSquared = (int) Math.round(Math.pow(radius*ratio, 2D));
  1406.         //int radiusSquared = radius*radius;
  1407.  
  1408.         int y1 = (int)Math.floor (y - radius);
  1409.         int y2 = (int)Math.ceil (y + radius) + 1;
  1410.         int x1 = (int)Math.floor (x - radius);
  1411.         int x2 = (int)Math.ceil (x + radius) + 1;
  1412.  
  1413.         if(y1 < 0)
  1414.             y1 = 0;
  1415.         if(y2 > height)
  1416.             y2 = height;
  1417.         if(x1 < 0)
  1418.             x1 = 0;
  1419.         if(x2 > width)
  1420.             x2 = width;
  1421.  
  1422.         for(int indexR = 0; indexR<depth; indexR++){
  1423.             for(int i = y1; i < y2; i++) {
  1424.                 //Check for interrupt
  1425.                 if(cancelThread) return;
  1426.                 for(int j = x1; j < x2; j++) {               
  1427.                     if((int) (Math.pow (j - x, 2D) + Math.pow (i - y, 2D)) < radiusSquared) {
  1428.                         houghValues[j][i][indexR] = 0;
  1429.                     }
  1430.                 }
  1431.             }
  1432.         }
  1433.     }
  1434.    
  1435.     /** Create a Thread[] array as large as the number of processors available.
  1436.     * From Stephan Preibisch's Multithreading.java class. See:
  1437.     * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
  1438.     */  
  1439.     private Thread[] newThreadArray() {  
  1440.         int n_cpus = Runtime.getRuntime().availableProcessors();  
  1441.         return new Thread[n_cpus];  
  1442.     }
  1443.    
  1444.     /** Start all given threads and wait on each of them until all are done.
  1445.     * From Stephan Preibisch's Multithreading.java class. See:
  1446.     * http://repo.or.cz/w/trakem2.git?a=blob;f=mpi/fruitfly/general/MultiThreading.java;hb=HEAD
  1447.      * @param threads
  1448.     */  
  1449.     public static void startAndJoin(Thread[] threads)  
  1450.     {  
  1451.         for (int ithread = 0; ithread < threads.length; ++ithread)  
  1452.         {  
  1453.             threads[ithread].setPriority(Thread.NORM_PRIORITY);  
  1454.             threads[ithread].start();  
  1455.         }  
  1456.  
  1457.         try  
  1458.         {    
  1459.             for (int ithread = 0; ithread < threads.length; ++ithread)  
  1460.                 threads[ithread].join();  
  1461.         } catch (InterruptedException ie)  
  1462.         {  
  1463.             throw new RuntimeException(ie);  
  1464.         }  
  1465.     }
  1466.    
  1467.     //Catches publish info from background thread so that the status can be checked
  1468.     @Override
  1469.     protected void process(List<String> status) {
  1470.       currentStatus = status.get(status.size()-1);
  1471.    }
  1472.    
  1473.    //Allows GUI class to get status updates from the worker thread
  1474.    public String getStatus(){
  1475.        return currentStatus;
  1476.    }
  1477.    
  1478.    //Flags that all active child threads should stop when cancel button is pressed in GUI class
  1479.    public void interruptThreads(boolean a){
  1480.        this.cancelThread = a;
  1481.    }
  1482. }