automatically draw a ROI circle

[Adam Hacking]

Hi,
 
I have an object of interest in an image that is round but has a discontinous perimeter (cross-section through a metallic foam rod). I would like to automatically draw a ROI (circle) that best fits the outline. In user mode this is simple as one simply constrains or shrinks the ROI until it touches the metal. The result is the cross section bound by a round ROI.
 
Can anyone suggest a simple approach to code this so it does not have to be done manually ?
 
Thanks
 
Adam

[Michael Schmid]

Hi Adam,

this question is not really clear to me.
- Inscribed circle: get it from Process>Binary>Ultimate Points
- Best-fit circle: get the centroid and area from the particle analyzer
- Circumscribed circle: modify the macro
   http://rsb.info.nih.gov/ij/macros/FeretsDiameter.txt
to determine the circumscribed circle and set the roi.

(BTW, the center of the circumscribed circle would be a nice addition  
to the Particle analyzer output...).

Michael
________________________________________________________________

On 10 Jun 2009, at 17:41, Adam Hacking wrote:

> Hi,
>
> I have an object of interest in an image that is round but has a  
> discontinous perimeter (cross-section through a metallic foam rod).  
> I would like to automatically draw a ROI (circle) that best fits  
> the outline. In user mode this is simple as one simply constrains  
> or shrinks the ROI until it touches the metal. The result is the  
> cross section bound by a round ROI.
>
> Can anyone suggest a simple approach to code this so it does not  
> have to be done manually ?
>
> Thanks
>
> Adam

[Gabriel Landini]

On Wednesday 10 June 2009 17:16:32 Michael Schmid wrote:
> this question is not really clear to me.
> - Inscribed circle: get it from Process>Binary>Ultimate Points
> - Best-fit circle: get the centroid and area from the particle analyzer
> - Circumscribed circle: modify the macro
>    http://rsb.info.nih.gov/ij/macros/FeretsDiameter.txt
> to determine the circumscribed circle and set the roi.

Or you can get the Minimum Bounding Circle with the Convex_Hull_Plus plugin.
http://www.dentistry.bham.ac.uk/landinig/software/software.html

> (BTW, the center of the circumscribed circle would be a nice addition
> to the Particle analyzer output...).

You can get the centre coordinates and the radius with the Particles8 plugin
(use the 'morphology' option).

Cheers

Gabriel

[Michael Schmid]

Below is a macro that draws the minimum enclosing circle for a selection
as, e.g. generated by the wand.
In case of problems, let me know.

Mind possible line breaks introduced by the mailer: Only the 'macro',
'var' and 'function' statements as well as '' and comments start in column
1, i.e., without preceding blanks. If there are other statements in column
1, you have to join these line with the previous ones.

Michael
___________________________________________________________________________

/* This macro creates a circular selection that is the smallest circle
   enclosing the current selection.
   Version: 2009-06-12 Michael Schmid

   Restrictions:
   - Does not work with composite selections
   - Due to rounding errors, some selection points may be slightly outside
the circle
*/


//global variables used for passing or as return values
var fourIndices = newArray(4);
var xcenter, ycenter, radius;

macro 'Smallest Enclosing Circle' {
  requires("1.29n");
  if (selectionType<0) exit("Error: Roi Required");
  if (selectionType==9) exit("This macro does not work with composite
selections");
  run("Line Width...", "line=1");
  getSelectionCoordinates(xCoordinates, yCoordinates);
  smallestEnclosingCircle(xCoordinates, yCoordinates);
  diameter = round(2*radius);
  makeOval(round(xcenter-radius), round(ycenter-radius), diameter, diameter);
}

/* Finds the smallest circle enclosing a set of points */
/* Input: arrays of x and y coordinates of the points */
/* Returns global variables xcenter, ycenter, radius */
function smallestEnclosingCircle(x,y) {
  n = x.length;
  if (n==1)
    return newArray(x[0], y[0], 0);
  else if (n==2)
    return circle2(x[0], y[0], x[1], y[1]);
  else if (n==3)
    return circle3(x[0], y[0], x[1], y[1], x[2], y[2]);
  //As starting point, find indices of min & max x & y
  xmin = 999999999; ymin=999999999; xmax=-1; ymax=-1;
  for (i=0; i<n; i++) {
    if (x[i]<xmin) {xmin=x[i]; fourIndices[0]=i;}
    if (x[i]>xmax) {xmax=x[i]; fourIndices[1]=i;}
    if (y[i]<ymin) {ymin=y[i]; fourIndices[2]=i;}
    if (y[i]>ymax) {ymax=y[i]; fourIndices[3]=i;}
  }
  do {
    badIndex = circle4(x, y);  //get circle through points listed in
fourIndices
    newIndex = -1;
    largestRadius = -1;
    for (i=0; i<n; i++) {      //get point most distant from center of circle
      r = vecLength(xcenter-x[i], ycenter-y[i]);
      if (r > largestRadius) {
        largestRadius = r;
        newIndex = i;
      }
    }
    //print(largestRadius);
    retry = (largestRadius > radius*1.0000000000001);
    fourIndices[badIndex] = newIndex; //add most distant point
  } while (retry);
}

//circle spanned by diameter between two points.
function circle2(xa,ya,xb,yb) {
  xcenter = 0.5*(xa+xb);
  ycenter = 0.5*(ya+yb);
  radius = 0.5*vecLength(xa-xb, ya-yb);
  return;
}
//smallest circle enclosing 3 points.
function circle3(xa,ya,xb,yb,xc,yc) {
  xab = xb-xa; yab = yb-ya; c = vecLength(xab, yab);
  xac = xc-xa; yac = yc-ya; b = vecLength(xac, yac);
  xbc = xc-xb; ybc = yc-yb; a = vecLength(xbc, ybc);
  if (b==0 || c==0 || a*a>=b*b+c*c) return circle2(xb,yb,xc,yc);
  if (b*b>=a*a+c*c) return circle2(xa,ya,xc,yc);
  if (c*c>=a*a+b*b) return circle2(xa,ya,xb,yb);
  d = 2*(xab*yac - yab*xac);
  xcenter = xa + (yac*c*c-yab*b*b)/d;
  ycenter = ya + (xab*b*b-xac*c*c)/d;
  radius = vecLength(xa-xcenter, ya-ycenter);
  return;
}
//Get enclosing circle for 4 points of the x, y array and return which
//of the 4 points we may eliminate
//Point indices of the 4 points are in global array fourIndices
function circle4(x, y) {
  rxy = newArray(12); //0...3 is r, 4...7 is x, 8..11 is y
  circle3(x[fourIndices[1]], y[fourIndices[1]], x[fourIndices[2]],
y[fourIndices[2]], x[fourIndices[3]], y[fourIndices[3]]);
  rxy[0] = radius; rxy[4] = xcenter; rxy[8] = ycenter;
  circle3(x[fourIndices[0]], y[fourIndices[0]], x[fourIndices[2]],
y[fourIndices[2]], x[fourIndices[3]], y[fourIndices[3]]);
  rxy[1] = radius; rxy[5] = xcenter; rxy[9] = ycenter;
  circle3(x[fourIndices[0]], y[fourIndices[0]], x[fourIndices[1]],
y[fourIndices[1]], x[fourIndices[3]], y[fourIndices[3]]);
  rxy[2] = radius; rxy[6] = xcenter; rxy[10] = ycenter;
  circle3(x[fourIndices[0]], y[fourIndices[0]], x[fourIndices[1]],
y[fourIndices[1]], x[fourIndices[2]], y[fourIndices[2]]);
  rxy[3] = radius; rxy[7] = xcenter; rxy[11] = ycenter;
  radius = 0;
  for (i=0; i<4; i++)
    if (rxy[i]>radius) {
      badIndex = i;
      radius = rxy[badIndex];
    }
  xcenter = rxy[badIndex + 4]; ycenter = rxy[badIndex + 8];
  return badIndex;
}

function vecLength(dx, dy) {
  return sqrt(dx*dx+dy*dy);
}