001 import ij.IJ;
002 import ij.ImagePlus;
003 import ij.WindowManager;
004 import ij.gui.OvalRoi;
005 import ij.gui.PolygonRoi;
006 import ij.gui.Roi;
007 import ij.measure.Measurements;
008 import ij.measure.ResultsTable;
009 import ij.plugin.filter.Analyzer;
010 import ij.plugin.filter.PlugInFilter;
011 import ij.process.ImageProcessor;
012 import ij.process.ImageStatistics;
013 import java.awt.Point;
014 import java.awt.Rectangle;
015 import java.awt.geom.Point2D;
016 /**
017 * Plugin for ImageJ.
018 * Finds and draws a convex hull
019 * by Thomas R. Roy, University of Alberta, Canada.
020 * Finds the smallest bounding circle using three points of
021 * the convex hull for a binary image, as proposed by Xavier Draye, and
022 * a modification in the case that this algorithm does not work, using
023 * the largest triangle made from points on the hull.
024 * The image has to be binary.
025
026 * @since date: 2002-2004
027 * @since updated January 1, 2005
028 * @since updated to check for binary images May 23, 2005
029 * @author A. Karperien
030 * Charles Sturt University, Australia
031 * @author Thomas R. Roy
032 * University of Alberta, Canada
033 * @author Wayne Rasband, NIH
034 */
035
036 public class HullAndCircle05_ implements PlugInFilter, Measurements {
037 /**variables for circularity*/
038 float majorcircle=0, minorcircle=0, midx=0, midy=0;
039 /**Variable for bounding circle's diameter.*/
040 double diameter=0;
041 private int width,height; //width and height of window
042 private int imgHeight, imgWidth; //dimensions of pixelated area
043 //private float centrex, centrey;
044 private int left, top, right, bottom;
045 private int topxc, rightyc, leftyc, bottomxc;
046 private String filename;
047 private ImagePlus img;
048 ImageStatistics stats=null;
049 int foreground = 0;
050 int X[] = new int [maxImageWidth * 2];
051 int Y[] = new int [maxImageWidth * 2];
052
053 private static final int maxImageWidth = 5000;
054 /**
055 * Returns the array of ints based on the doubles.
056 *
057 * @param A array of doubles
058 * @param length length of array to convert
059 * @return int array
060 */
061 public static int[] DoubleArrayToIntArray(double[] A, int length) {
062 int[]N=new int[length];
063 for (int i = 0; i < length; i++)
064 N[i]=(int)(A[i]);
065 return N;
066 }
067 /**
068 * Gets the centre and diameter for
069 * the circle from the passed three points
070 * (a,b), (c,d), and (e,f). Solves three equations
071 * simultaneously using the basic
072 * equation for a circle.
073 * (h-x)<sup>2</sup>+(k-y)<sup>2</sup>=r<sup>2</sup>
074 * Does not check for colinearity.
075 * @param a x coordinate of first point
076 * @param b y coordinate of first point
077 * @param c x coordinate of second point
078 * @param d y coordinate of second point
079 * @param e x coordinate of third point
080 * @param f y coordinate of third point
081 * @return double [] containing
082 * the centre coordinates(h, k) and the diameter
083 * of the circle
084 */
085 public double[] getCircle(double a, double b,
086 double c, double d, double e, double f) {
087 b=-b;
088 d=-d;
089 f=-f;
090
091 double k = ((sq(a)+sq(b))*(e-c)
092 + (sq(c)+sq(d))*(a-e) + (sq(e)+sq(f))*(c-a)) /
093 (2*(b*(e-c)+d*(a-e)+f*(c-a))) ;
094 double h = ((sq(a)+sq(b))*(f-d) + (sq(c)+sq(d))*(b-f) +
095 (sq(e)+sq(f))*(d-b)) / (2*(a*(f-d)+c*(b-f)+e*(d-b)));
096
097 double r = Math.sqrt(sq(a-h) + sq(b-k));
098
099 return new double []{h, -k, 2*r};
100
101 } /**
102 * Returns the bounding circle for the binary image.
103 *
104 * Finds three points on the passed convex hull then
105 * calls the {@link #getCircle getCircle} method
106 * @param XX int array of x coordinates for the convex hull
107 * @param YY int array of y coordinates for the convex hull
108 * @param ip ImageProcessor of binary image that the convex
109 * hull is for
110 * @return OvalRoi the bounding circle
111 */
112 public OvalRoi FindBoundingCircle(int[] XX, int[] YY, ImageProcessor ip) {
113
114 //find the distance between the first two points
115 double max = Math.abs(Point2D.distance(
116 (double)XX[0], (double)YY[0], (double)XX[0], (double)YY[0]));
117
118 double oldmax=max;
119
120 int max_i = 0;
121 int max_i2= 0;
122 double min = 0;
123 //find the greatest distance between two points on the
124 //convex hull
125 for (int anchor=0; anchor<XX.length; anchor++) {
126 for (int i = 0; i< XX.length; i++) {
127 double dist = Math.abs(Point2D.distance((double)XX[i],
128 (double)YY[i], (double)XX[anchor], (double)YY[anchor]));
129 max = Math.max(dist, max);
130 if (max>oldmax)
131 {oldmax=max;
132 max_i=i;max_i2=anchor;}
133
134 }}
135 //draws the line
136 double x1=XX[max_i];
137 double y1=YY[max_i];
138 double x2=XX[max_i2];
139 double y2=YY[max_i2];
140 // ip.drawLine((int)x1, (int)y1, (int)x2, (int)y2);
141 majorcircle=(float)max;
142 //finds the midpoint
143 midx = (float)((x1+x2)/2.0);
144 midy =(float)((y1+y2)/2.0);
145 // ip.drawDot((int)midx, (int)midy);
146
147 //finds the maximum distance from the midpoint of
148 //that line to another point in the convex hull
149 max = Math.abs(Point2D.distance((double)XX[0],
150 (double)YY[0], midx, midy));
151 oldmax=max;
152 min=max;
153 int thirdpoint = max_i;
154 for (int i = 0; i< XX.length; i++) {
155
156 {double dist = Math.abs(Point2D.distance((double)XX[i],
157 (double)YY[i], midx, midy));
158 max = Math.max(dist, max);
159 min = Math.min(dist, min);
160 if (max>oldmax)
161 {oldmax=max;
162 thirdpoint=i;
163 }
164 }
165 }
166 double thirdx=XX[thirdpoint];
167 double thirdy=YY[thirdpoint];
168 boolean isdiameter=false;
169 if((max_i==thirdpoint)||(thirdpoint==max_i2)) {
170 isdiameter=true;
171 }
172 //decides if the point is on the current line
173
174 //draws a line from the midpoint
175 //of the first line to a next furthest point on the convex hull
176 //ip.drawLine((int)midx, (int)midy, (int)XX[thirdpoint],
177 //(int)YY[thirdpoint]);
178 minorcircle=(float)max;
179 //Gets the bounding circle from the three points
180 //just found then draws it
181 double [] circ = getCircle(XX[max_i], YY[max_i],
182 XX[max_i2], YY[max_i2], thirdx, thirdy);
183 int [] Circ =
184 DoubleArrayToIntArray(circ, circ.length);
185 //IJ.showMessage("A"+Circ[0]+"A"+Circ[1]+"A"+Circ[2]);
186 OvalRoi boundingCircle = null;
187
188 if (!isdiameter) {
189 diameter=Circ[2];
190 boundingCircle=
191 new OvalRoi(Circ[0]-Circ[2]/2,
192 Circ[1]-Circ[2]/2,
193 Circ[2], Circ[2]);
194 midx=Circ[0]-Circ[2]/2;
195 midy=Circ[1]-Circ[2]/2;
196 }
197 else {
198 diameter=(float)(2.0*max);
199 boundingCircle=
200 new OvalRoi((int)midx-(int)(max), (int)midy-(int)(max),
201 (int)(2*max), (int)(2*max));
202 midx = (int)midx-(int)(max);
203 midy=(int)midy-(int)(max);
204 }
205
206 //ip.drawPolygon(boundingCircle.getPolygon());
207 return boundingCircle;
208 }
209
210
211
212 /**
213 * ImageJ set up
214 * @param arg arguments to run plugin = null
215 * @param imp ImageProcessor on which to work
216 * @return int an integer
217 */
218 public int setup(String arg, ImagePlus imp) {
219 if (imp==null)
220 return(DONE);
221 else { stats=imp.getStatistics();
222 IJ.register(this.getClass());
223 filename=imp.getTitle();
224 return IJ.setupDialog(imp, DOES_8G);
225 }
226 }
227
228
229
230 /**
231 * Draws the convex hull for the current ImageProcessor,
232 * makes some measurements, and draws the bounding circle
233 * on the image. The bounding circle is currently
234 * drawn only; to access it, use the BoundingCircle
235 * roi made in this method.
236 * @param ip ImageProcessor on which to work
237 * @see GetDimensions
238 * @see getPolygon
239 * @see FindTriangle
240 */
241
242 public void run(ImageProcessor ip) {
243 if (true) {
244 if (stats.histogram[0]>stats.histogram[255])
245 foreground = 255;
246 else foreground = 0;
247 {//ImageStatistics stats = imp.getStatistics();
248 if (stats.histogram[0]+stats.histogram[255]!=stats.pixelCount) {
249 IJ.error("8-bit binary image (0 and 255) required.");
250 return;
251 }
252
253 }
254 GetDimensions(ip);
255 PolygonRoi poly = this.getPolygon(ip);
256 img.setRoi(poly);
257 int ms = Analyzer.getMeasurements();
258 ms |= AREA+PERIMETER+CENTER_OF_MASS;//;+ELLIPSE;
259 Analyzer.setMeasurements(ms);
260 Analyzer ana = new Analyzer();
261 stats = img.getStatistics(ms);
262 Roi roi = img.getRoi();
263 ana.saveResults(stats, roi);
264 ResultsTable rt =Analyzer.getResultsTable();
265
266 ip.drawPolygon(poly.getPolygon());
267 OvalRoi BoundingCircle = FindTriangle(X, Y, ip, poly);
268 ip.drawPolygon(BoundingCircle.getPolygon());
269
270 int cr = rt.getCounter();
271 double area = rt.getValue(ResultsTable.AREA, cr-1);
272 double perimeter = rt.getValue(ResultsTable.PERIMETER, cr-1);
273 double circul = (perimeter*perimeter)/area;
274 rt.addValue("Circularity",
275 perimeter==0.0?0.0:4.0*Math.PI*(area/(perimeter*perimeter)));
276 rt.addValue("Roundness", circul);
277 rt.addLabel("FileName", filename);
278 // float major = rt.getValue(ResultsTable.MAJOR, cr-1);
279 // float minor = rt.getValue(ResultsTable.MINOR, cr-1);
280 // rt.addValue("Axis1", major);
281 // rt.addValue("Axis2", minor);
282 img.killRoi();
283 //Rectangle rect =poly.getBoundingRect();
284 //// ip.drawPolygon(poly.getPolygon());
285 //// OvalRoi BoundingCircle = FindTriangle(X, Y, ip, poly);
286 //// ip.drawPolygon(BoundingCircle.getPolygon());
287 rt.addValue("Bounding Circle Diameter", diameter);
288 ana.displayResults();
289 ana.updateHeadings();
290 }
291 }
292
293 /**Returns the square of the passed variable.
294 *@param N a double
295 *@return N*N
296 */
297 public double sq(double N){ return (N*N);}
298
299
300
301
302
303
304
305
306 /**
307 * This method determines the vertices of a polygon which forms a convex
308 * hull around the image. It does this by first mapping the vertical
309 * contours of the image. (Horizontal or vertical contours could have
310 * been used, but only one is necessary.) After they are mapped, all of
311 * the points which can be enclosed by connecting any other points are
312 * eliminated. This leaves only the outermost points, the ones that form
313 * the convex hull. Returns a polygon created with these points.
314 * @param ip ImageProcessor to work on
315 * @return PolygonRoi the convex hull
316 *@see mapTopContours
317 *@see mapBottomContours
318 *@see removeContours
319 *@author Thomas R. Roy, University of Alberta
320 */
321 public PolygonRoi getPolygon(ImageProcessor ip){
322
323 Point coords[] = new Point[maxImageWidth * 2];
324 for (int i = 0; i < maxImageWidth * 2; i++)
325 coords[i] = new Point();
326 int[] histogram = new int[256];
327 int vertCount=0;
328 int y, x;
329 Point farLeft, farRight;
330
331 //Start by mapping farthest left point
332 x = this.left;
333 y=top-2;
334 do{
335 y++;
336 ip.setRoi(x, y, 1, 1);
337 histogram = ip.getHistogram();
338 } while (( histogram[ this.foreground ]==0 ) && y<=bottom);
339 farLeft = new Point(x,y);
340 int index = x;
341
342 coords[0].x=index;
343 coords[0].y=y;
344
345 //Find the furthest right point
346 x = this.right;
347 y=top-2;
348 do{
349 y++;
350 ip.setRoi(x, y, 1, 1);
351 histogram = ip.getHistogram();
352 } while (( histogram[ this.foreground ]==0 ) && y<=bottom);
353 farRight = new Point(x,y);
354
355 //Map the top and bottom contours
356 vertCount = mapTopContours(farLeft, farRight, coords, vertCount, ip);
357 int endPoint = vertCount;
358 vertCount = mapBottomContours(farLeft, farRight, coords, vertCount, ip);
359
360 vertCount++;
361 coords[vertCount]=coords[0];
362
363 //Remove the concave contours, leaving only concave ones
364 coords = removeContours(coords, vertCount, endPoint);
365
366 //Produce an array for the output
367 X = new int [coords.length];
368 Y = new int [coords.length];
369 for (int i = 0; i < coords.length; i++) {
370 X[i]=coords[i].x;
371 Y[i]=coords[i].y;
372 }
373 return new PolygonRoi(X, Y, coords.length, img,Roi.POLYGON);
374 }
375 /**Maps the bottom contours of the image.
376 @param farLeft Point
377 @param farRight Point
378 @param coords Point[]
379 @param vertCount int
380 @param ip ImageProcessor
381 @return int for vertical count*/
382 public int mapBottomContours(Point farLeft, Point farRight,
383 Point[] coords, int vertCount, ImageProcessor ip){
384 //Begin mapping bottom contours, right to left
385 int[] histogram = new int[256];
386
387 for (int index = farRight.x; index >= farLeft.x; index--){
388 int y = this.bottom;
389 do{
390 y--;
391 ip.setRoi(index, y, 1, 1);
392 histogram = ip.getHistogram();
393 } while (( histogram[ this.foreground ]==0 ) && y>=top);
394
395 if (y>top){
396 vertCount++;
397 coords[vertCount].x = index;
398 coords[vertCount].y=y;
399 }
400 }
401 return vertCount;
402
403 }
404 /**Maps top contours.
405 *@param farLeft Point
406 *@param farRight Point
407 *@param coords Point[]
408 *@param vertCount int
409 *@param ip ImageProcessor
410 *@return int for top count*/
411 public int mapTopContours(Point farLeft, Point farRight,
412 Point[] coords, int vertCount, ImageProcessor ip){
413 //Begin mapping contours left to right
414 int[] histogram = new int[256];
415
416 for (int index = farLeft.x; index <= farRight.x; index++){
417
418 int y = this.top;
419 do{
420 y++;
421 ip.setRoi(index, y, 1, 1);
422 histogram = ip.getHistogram();
423 } while (( histogram[ this.foreground ]==0 ) && y<=(bottom));
424
425 if (y<bottom){
426 vertCount++;
427 coords[vertCount].x = index;
428 coords[vertCount].y = y;
429 }
430 }
431 return vertCount;
432
433 }
434 /**
435 * This method removes contours that contain concave areas.
436 * The algorithm tests each contour line to ensure that it has a greater
437 * slope than all of the lines after it. If this is not the case, then it
438 * replaces the endpoint with the endpoint of the next contour pixel, and
439 * checks again, until it has removed all of the contour points that are
440 * concave.
441 * @return Point[] array of new coordinates
442 */
443 public Point[] removeContours(Point coords[], int length, int furthest){
444
445 int index=0;
446 int nIndex = 0;
447 float rise, run, slope, slope2;
448
449 //First check the top contours
450 do{
451 float highest = 0;
452 int highestIndex = -1;
453 for(int i = index+1; i <= furthest; i++){
454 run = coords[i].x - coords[index].x;
455 rise = coords[i].y - coords[index].y;
456 if (run == 0)
457 slope = 100000;
458 else
459 slope = rise/run;
460
461 if (slope <= highest || highestIndex == -1){
462 highest = slope;
463 highestIndex = i;
464 }
465 }
466 coords[nIndex] = coords[index];
467 nIndex++;
468 index = highestIndex;
469
470 } while (index != -1);
471
472 //Now check the bottom contours
473 index=furthest;
474 do{
475 float highest = 0;
476 int highestIndex = -1;
477 for(int i = index+1; i < length; i++){
478 run = coords[i].x - coords[index].x;
479 rise = coords[i].y - coords[index].y;
480 if (run == 0)
481 slope = 100000;
482 else
483 slope = rise/run;
484
485 if (slope <= highest || highestIndex == -1){
486 highest = slope;
487 highestIndex = i;
488 }
489 }
490 coords[nIndex] = coords[index];
491 nIndex++;
492 index = highestIndex;
493
494 } while (index != -1);
495
496 Point newCoords[] = new Point[nIndex];
497 for (int i = 0; i < nIndex; i++)
498 newCoords[i] = coords[i];
499
500 return newCoords;
501 }
502
503
504 /**Determine the dimensions of the image,
505 * and sets some of the important variables.
506 *@param ip ImageProcessor on which to operate
507 */
508 public void GetDimensions(ImageProcessor ip) {
509
510
511 img = WindowManager.getCurrentImage();
512 int[] histogram = new int[256];/*make an array for histogram*/
513 this.width = ip.getWidth();/*record the width and height*/
514 this.height = ip.getHeight();
515 Rectangle roi;
516
517 //start at the leftmost edge, move to the middle,
518 //and find the first coloured pixel for the left edge
519 this.left = -1;
520 do
521 { this.left++;
522 ip.setRoi(this.left, 0, 1, this.height);
523 histogram = ip.getHistogram();
524 } while ( histogram[this.foreground]==0 );
525
526 //start at the top and find the first pixel
527 this.top = -1;
528 do
529 { this.top++;
530 ip.setRoi(this.left, this.top, this.width, 1);
531 histogram = ip.getHistogram();
532 } while (histogram[this.foreground]==0);
533
534 //start at the right and find the first pixel
535 right =width-1;
536 do
537 { this.right--;
538 ip.setRoi(this.right, this.top, 1, this.height);
539 histogram = ip.getHistogram();
540 } while (histogram[this.foreground]==0);
541
542 //start at the bottom of the image and find the first pixel*/
543 this.bottom =this.height;
544 do
545 { this.bottom--;
546 ip.setRoi(this.left, this.bottom, this.right, 1);
547 histogram = ip.getHistogram();
548 } while (histogram[this.foreground]==0);
549 // This part of the methods finds the coordinates
550 //of the top, bottom, right and left pixels,
551 // which is used later for the polygon.
552 this.leftyc = this.top-1;
553 do
554 { this.leftyc++;
555 ip.setRoi(this.left, this.leftyc, 1, 1);
556 histogram = ip.getHistogram();
557 } while (histogram[this.foreground]==0);
558
559 this.topxc = this.left-1;
560 do
561 { topxc++;
562 ip.setRoi(topxc, top, 1, 1);
563 histogram = ip.getHistogram();
564 } while (histogram[this.foreground]==0);
565
566 this.rightyc = this.top-1;
567 do
568 { this.rightyc++;
569 ip.setRoi(this.right, this.rightyc, 1, 1);
570 histogram = ip.getHistogram();
571 } while (histogram[this.foreground]==0);
572
573 this.bottomxc = this.right+1;
574 do
575 { this.bottomxc--;
576 ip.setRoi(this.bottomxc, this.bottom, 1, 1);
577 histogram = ip.getHistogram();
578 } while (histogram[this.foreground]==0);
579
580 //We have to ensure these regions don't overlap.
581 this.top--;
582 this.bottom++;
583
584 this.imgHeight = bottom - top;
585 this.imgWidth = right - left;
586 //this.centrex = (float)(left + (this.imgWidth/2));
587 //this.centrex = (float)(top + (this.imgHeight/2));
588 }
589
590 /**
591 * Finds the triangle having the greatest perimeter using
592 * all points in the convex hull. Then uses this as the three
593 * points on the hull that define the circle. Is an alternative
594 * when Draye's method yields a triangle that does not include all
595 * of the points.
596 * @return OvalRoi a circle as an OvalRoi
597 * @param XX array of x coordinates for the hull
598 * @param YY array of y coordinates for the hull
599 * @param ip ImageProcessor on which to work
600 * @param poly PolygonRoi for the convex hull
601 * @see FindBoundingCircle
602 */
603 public OvalRoi FindTriangle(int[] XX, int[] YY,
604 ImageProcessor ip, PolygonRoi poly) {
605 //find perimeter of each triangle
606
607 OvalRoi boundingCircle = null;
608 boundingCircle=FindBoundingCircle(XX, YY, ip);
609 boolean AllIn=true;
610 for (int i = 0; i < XX.length; i++) {
611 if
612 (!(boundingCircle.getPolygon().contains(new
613 Point(poly.getXCoordinates()[i],
614 poly.getYCoordinates()[i]))))
615 if(!boundingCircle.getPolygon().intersects(XX[i]-1, YY[i]-1,
616 3,3))
617 {AllIn=false;
618 //IJ.showMessage("not at"+i);
619 }
620 }
621 if (AllIn==false)
622 { double p = 0;
623 //find the distance between the first two points
624 double oldmin = Math.abs(Point2D.distance(
625 (double)XX[0], (double)YY[0], (double)XX[0], (double)YY[0]));
626
627 double min=oldmin;
628
629 int max_i = 0;
630 int max_i2= 0;
631
632 //find the greatest distance between two points on the
633 //convex hull
634 for (int anchor=0; anchor<XX.length; anchor++) {
635 for (int i = 0; i< XX.length; i++) {
636 double dist = Math.abs(Point2D.distance((double)XX[i],
637 (double)YY[i], (double)XX[anchor], (double)YY[anchor]));
638 min = Math.max(dist, min);
639 if (min>oldmin)
640 {oldmin=min;
641 max_i=i;max_i2=anchor;}
642
643 }}
644 int I=0, J=0, A=0;
645 double maxdist=0;
646
647 for (int anchor=0; anchor<XX.length; anchor++) {
648 for (int i = 0; i< XX.length; i++) {
649 for (int j = 0; j< XX.length; j++) {
650 double dist1 = Math.abs(Point2D.distance((double)XX[i],
651 (double)YY[i], (double)XX[anchor], (double)YY[anchor]));
652 double dist2 = Math.abs(Point2D.distance((double)XX[j],
653 (double)YY[j], (double)XX[anchor], (double)YY[anchor]));
654 double dist3 = Math.abs(Point2D.distance((double)XX[i],
655 (double)YY[i], (double)XX[j], (double)YY[j]));
656 if (true)//((dist1>=min)||(dist2>=min)||(dist3>=min))
657 { double newdist=dist1+dist2+dist3;
658 if (newdist>maxdist){
659 maxdist=newdist;
660 I=i;
661 J=j;
662 A=anchor;
663 double [] circ = getCircle(XX[I], YY[I],
664 XX[J], YY[J], XX[A], YY[A]);
665 int [] Circ =
666 DoubleArrayToIntArray(circ, circ.length);
667 //IJ.showMessage("A"+Circ[0]+"A"+Circ[1]+"A"+Circ[2]);
668
669
670 if (true) {
671 diameter=Circ[2];
672 boundingCircle=
673 new OvalRoi(Circ[0]-Circ[2]/2,
674 Circ[1]-Circ[2]/2,
675 Circ[2], Circ[2]);
676 midx=Circ[0]-Circ[2]/2;
677 midy=Circ[1]-Circ[2]/2;
678 }
679 } }
680 }
681 }
682 }
683
684 }
685 return boundingCircle;
686 }
687
688
689
690 }//end class