Fit ellipse in the middle of two ellipses

Hi everyone,

I have two ellipses that I drew in ImageJ (a bigger one around a smaller one). Now I want to automatically place an ellipse exactly in between, so that the distance of the middle ellipse to the outer and inner ellipse is equal. Does anyone know of a macro or function to do this?

Thanks!

Soufiane

1 Like

Hi @soufiane,

there is a function which can tell you the position and size of the current ROI called Roi.getBounds, you can explore what ROI methods are there by starting to type Roi. in the script editor.

Furthermore, there are make... functions for generating new ROIs:

I you need more help in programming your macro, a screenshot of your ellipses might be helpful to understand where the new ROI is supposed to be :wink:

Cheers,
Robert

Hi Robert,

Thank you for your response.
I would appreciate if you could help me with programming the macro.
I have the image below, to give an impression of what I would like.

Hi
@soufiane
Do your two ellipses have (- same center? - same directions of the two axes?)

Hi @Mathew,

I don’t think so. First I have drawn polygons and after that I have done ‘Fit Ellipse’.

Edit:
I now realize there should be a center to make an ellipse in between. I think we could use the center of the inner circle.

Hi
@soufiane
Is this okay?

You can adapt the macro

setBatchMode(true);
setOption("BlackBackground", true);
newImage("Untitled", "8-bit black", 512, 512, 0);
makePolygon(123,304,148,206,227,168,307,148,357,148,419,216,363,267,262,332,155,358,122,355,114,330);
run("Fit Ellipse");
roiManager("Add");
//setTool("polygon");
makePolygon(273,108,147,132,97,260,53,368,71,427,154,446,305,441,421,335,443,210,473,96,453,54,376,30,297,54);
run("Fit Ellipse");
roiManager("Add");
run("Set Measurements...", "area centroid add redirect=None decimal=0");
roiManager("Select", 0);
roiManager("Measure");
x=getResult("X",0);
y=getResult("Y",0);
makePoint(x,y);
roiManager("Add");
setBackgroundColor(255, 255, 255);
roiManager("Select", 1);
run("Clear Outside");
roiManager("Select", 0);
run("Clear", "slice");
run("Select None");
run("Convert to Mask");
run("Skeletonize");
roiManager("Show All without labels");
setBatchMode(false);
exit();
2 Likes

Hello Mathew (and Soufiane) -

The problem is that the result of Skeletonize need not be an
ellipse. So if you want an ellipse in between the two original
ellipses, this won’t work.

To Soufiane: Do you just need a smooth curve in between your
two ellipses, or does it have to be another ellipse?

Thanks, mm

Hi
@mountain_man

Beautiful observation !
I will then add a threshold on the skeleton + Fit ellipse.
It gives that.

1 Like

Hi @mountain_man
A smooth curve would be perfect indeed.
The one @Mathew showed in his last post is perfect. What did you add to the previous macro?

Thanks

Edit:
I only get the center point with the previous macro, not the wite curve in between.

Hi
@soufiane
I will then add a threshold on the skeleton + Fit ellipse.

setBatchMode(true);
setOption("BlackBackground", true);
newImage("Untitled", "8-bit black", 512, 512, 0);
run("Duplicate...", " ");
makePolygon(123,304,148,206,227,168,307,148,357,148,419,216,363,267,262,332,155,358,122,355,114,330);
run("Fit Ellipse");
roiManager("Add");
//setTool("polygon");
makePolygon(273,108,147,132,97,260,53,368,71,427,154,446,305,441,421,335,443,210,473,96,453,54,376,30,297,54);
run("Fit Ellipse");
roiManager("Add");
run("Set Measurements...", "area centroid add redirect=None decimal=0");
roiManager("Select", 0);
roiManager("Measure");
x=getResult("X",0);
y=getResult("Y",0);
makePoint(x,y);
roiManager("Add");
setBackgroundColor(255, 255, 255);
roiManager("Select", 1);
run("Clear Outside");
roiManager("Select", 0);
run("Clear", "slice");
run("Select None");
run("Convert to Mask");
run("Skeletonize");
setAutoThreshold("Default dark");
//run("Threshold...");
run("Convert to Mask");
run("Analyze Particles...", "display add in_situ");
roiManager("Select", 3);
run("Fit Ellipse");
roiManager("Add");
roiManager("Select", 3);
roiManager("Delete");
selectWindow("Untitled");
roiManager("Show All without labels");
roiManager("Select", 3);
roiManager("Set Color", "green");
roiManager("Set Line Width", 0);
run("Select None");
setBatchMode(false);
exit();
1 Like

Hi @Mathew, It runs perfect when I run it like it is.
Now I have written a short piece of macro where I would like to integrate your macro.
I use an RGB image and ask the user to draw two polygons which are fitted into ellipses and then I would like your macro to create the ellipse in between. I don’t get it to work, I don’t know why.
This is my macro

//Ask for user to place multi point for outer circle
waitForUser(“Draw polygon for outer circle”);
getSelectionCoordinates(x,y)
run(“Fit Ellipse”);
roiManager(“Add”);
roiManager(“Select”, 0);
roiManager(“Rename”, “Outer circle”);
roiManager(“Show All”);
roiManager(“Show None”);
roiManager(“Show All”);

//Ask for user to place multi point for inner circle
waitForUser(“Draw polygon for inner circle”);
getSelectionCoordinates(x,y)
run(“Fit Ellipse”);
roiManager(“Add”);
roiManager(“Select”, 1);
roiManager(“Rename”, “Inner circle”);
roiManager(“Show All”);
roiManager(“Show None”);
roiManager(“Show All”);

1 Like

I would define the largest ROI, then use “Edit>Selection>Enlarge” and set a negative value to reduce the ROI in size (or viceversa, with positive values to enlarge it).

2 Likes

Thanks @gabriel
I prefer the method of Mathew, since that prevents any problems when the two ellipses are differently orientated like in his example image.

Hello Mathew (and Soufiane) -

The problem is that the ellipse obtained by running Fit Ellipse
on the Skeletonize curve that is between the two original ellipses
need not itself be between those two ellipses.

I modified your macro to change the first polygon and therefore
the smaller of the two ellipses. This macro illustrates a case
where the fit ellipse does not entirely enclose the “inner” ellipse:

setBatchMode(true);
setOption("BlackBackground", true);
newImage("Untitled", "8-bit black", 512, 512, 0);
run("Duplicate...", " ");
// makePolygon(123,304,148,206,227,168,307,148,357,148,419,216,363,267,262,332,155,358,122,355,114,330);
makePolygon (110,270, 125,240, 202,265, 282,278, 332,280, 385,345, 338,360, 237,360, 130,330, 105,325, 100,300);
run("Fit Ellipse");
roiManager("Add");
//setTool("polygon");
makePolygon(273,108,147,132,97,260,53,368,71,427,154,446,305,441,421,335,443,210,473,96,453,54,376,30,297,54);
run("Fit Ellipse");
roiManager("Add");
run("Set Measurements...", "area centroid add redirect=None decimal=0");
roiManager("Select", 0);
roiManager("Measure");
x=getResult("X",0);
y=getResult("Y",0);
makePoint(x,y);
roiManager("Add");
setBackgroundColor(255, 255, 255);
roiManager("Select", 1);
run("Clear Outside");
roiManager("Select", 0);
run("Clear", "slice");
run("Select None");
run("Convert to Mask");
run("Skeletonize");
setAutoThreshold("Default dark");
//run("Threshold...");
run("Convert to Mask");
run("Analyze Particles...", "display add in_situ");
roiManager("Select", 3);
run("Fit Ellipse");
roiManager("Add");
roiManager("Select", 3);
roiManager("Delete");
selectWindow("Untitled");
roiManager("Show All without labels");
roiManager("Select", 3);
roiManager("Set Color", "green");
roiManager("Set Line Width", 0);
run("Select None");
setBatchMode(false);
exit();

Here is the resulting image:

To Soufiane: Is your use case such that it would be okay for your
“middle” ellipse to overlap the others? Or is your use case such
that your inner and outer ellipses are likely to be well enough aligned
with one another that this overlap is unlikely to occur in practice?

Or, as I think you said above, is Mathew’s Skeletonize curve
adequate, even though it isn’t elliptical?

Thanks, mm

Thank you @mountain_man
Very illustrative.
In my case, I think I can do both. The skeletonize result would be perfectly fine with me and like you said, the overlap you showed is unlikely to occur in my set of images.

I understand that it is mathematically not possible to get a middle ellipse in you figure?

Thank you both for your help!

Hello Soufiane -

As an aside:

I think it is mathematically possible to get the “in-the-middle”
ellipse.

For example, if you have an inner and outer ellipse (that are
not touching), you could always enlarge the inner ellipse just
a little bit – not enough to bump into the outer ellipse. Or you
could shrink the outer ellipse a little. Both of these would be
“in between” the inner and outer ellipses. But neither would
really be “in the middle.”

I think what you want is the ellipse that is in some sense
halfway between the inner and outer ellipses.

My problem is that I don’t know how to define what it means
to be “halfway between” in precise terms. I believe that there
is some nice mathematical definition of this – probably even a
canonical definition (that is, one clearly best definition) – I just
don’t know what it is.

It’s easy enough to give definitions for curves that are halfway
between the two ellipses. (Mathew’s Skeletonize curve is a
reasonable example of this.) I just don’t see any way to modify
those definitions work in the case where the curve is required
to be an ellipse.

Thanks, mm

1 Like

thanks.I learned a lot.

Hello everyone!

This topic is really interesting!

I was caught by @mountain_man’s reflections about what “in the middle” really means in this context, or what mathematical form the resulting answer might take.

As I’m kind of a sucker for this kind of stuff, I couldn’t resist to try some options [disclaimer: this might drift a little from @soufiane’s original goal, but I think that in the end it will contribute to overall discussion].

My first approach was based on interpolation: get the parameters that define both ellipses (center, angle, major and minor axis) and interpolate between the two sets to define an intermediate ellipse, in the same way one would implement a “morph” from one to the other. This seemed to produce reasonable results when trying ellipse pairs similar to those presented by @soufiane in the first example image, but then it performed poorly when tested with extreme examples like the one presented in this post by @mountain_man (where the inner ellipse is not oriented in the same direction as the outer one):

It looked like there were no reasonable candidates of “middle” ellipses that could work in this situation.

Here’s the code to try it:

if(!isOpen("test ellipses")){
	newImage("test ellipses", "8-bit black", 512, 512, 1);
	if(roiManager("count")>0){
		roiManager("deselect");
		roiManager("delete");
		}
	makeEllipse(71, 380, 406, 52, 0.50);
	roiManager("Add");
	makeEllipse(346, 55, 353, 205, 0.45);
	roiManager("Add");
	}

roiManager("select", 0);
List.setMeasurements;
cx1 = List.getValue("X");
cy1 = List.getValue("Y");
minor1 = List.getValue("Minor");
major1 = List.getValue("Major");
angle1 = -List.getValue("Angle")/180*PI;

roiManager("select", 1);
List.setMeasurements;
cx2 = List.getValue("X");
cy2 = List.getValue("Y");
minor2 = List.getValue("Minor");
major2 = List.getValue("Major");
angle2 = -List.getValue("Angle")/180*PI;


major = 0.5 * (major2-major1) + major1;
minor = 0.5 * (minor2-minor1) + minor1;
angle = 0.5 * (angle2-angle1) + angle1;
cx = 0.5 * (cx2-cx1) + cx1;
cy = 0.5 * (cy2-cy1) + cy1;

roiManager("Show All without labels");
makeEllipse(cx-cos(angle)*major/2, cy-sin(angle)*major/2, cx+cos(angle)*major/2, cy+sin(angle)*major/2, minor/major);
run("Properties... ", "stroke=orange width=2 fill=none");

I then went back to the definition of “in the middle”, and looked for general solution agnostic of the expected shape: the resulting curve should be -at any point- equidistant to both ellipses. This can be accomplished by finding pairs of points, one from each ellipse, that are the closest to each other, and then finding the mid-point of their connecting segment. The solution should be the geometrical locus of all these “middle points”. When I tried this approach with the example shown before, the resulting solution was certainly not looking like an ellipse:



Here’s the second script to reproduce this:

if(!isOpen("test ellipses")){
	newImage("test ellipses", "8-bit black", 512, 512, 1);
	if(roiManager("count")>0){
		roiManager("deselect");
		roiManager("delete");
		}
	makeEllipse(71, 380, 406, 52, 0.50);
	roiManager("Add");
	makeEllipse(346, 55, 353, 205, 0.45);
	roiManager("Add");
	}

run("Select None");

roiManager("select", 0);
getSelectionCoordinates(xp1, yp1);


roiManager("select", 1);
run("Fit Spline");
run("Interpolate", "interval=1 smooth");
getSelectionCoordinates(xp2, yp2);


xp3=newArray(xp1.length);
yp3=newArray(xp1.length);

for (i=0; i<xp1.length; i++){
	jmin=-1;
	d2min=-1; 
	for(j=0; j<xp2.length; j++){
		d2 = pow( xp1[i]-xp2[j], 2) + pow (yp1[i]-yp2[j], 2);
		if(d2<=d2min || d2min==-1){ 
			jmin=j;
			d2min=d2;
			}
		}
	xp3[i] = (xp1[i] + xp2[jmin]) /2;
	yp3[i] = (yp1[i] + yp2[jmin]) /2;
	}

roiManager("Show All without labels");
makeSelection("polygon", xp3, yp3);

run("Properties... ", "stroke=orange width=2 fill=none");

In other words, it became evident that the “middle curve” between two ellipses, one enclosing the other, is not an ellipse for the general case. As a corollary, any ellipse-fitting approach will perform worse the more dissimilar both ellipses become.

As for the exact mathematical form of this equidistant solution, and observing the resulting shapes, it looks to me that it will be way more involved than I expected at first. Perhaps someone can provide a pointer to an analytical aproach to these kind of problems. It certainly escapes my current skills.

Anyway, those are my two cents.

Thanks everybody for the discussion!

Cheers!
Nico

2 Likes

To all …
Thank you for this discussion.

Hello Soufiane and Mathew and Nico -

Upon further reflection, there does seem to be a systematic
way to construct the “in-the-middle” ellipse.

Here are some illustrative images:

In each image, the two thicker red ROIs are the inner and
outer ellipses. The thicker blue ROI is the “in-the-middle”
ellipse, and the other ellipses fill out the series of ellipses
that interpolate from the inner to the outer.

I’ve chosen “in the middle” to mean halfway between
in “linear scale,” where linear scale is
sqrt (major-axis * minor-axis) (proportional
to sqrt (area)).

Here is some text output from the script that generated the
images:

1.52t99
ellipses A:
linear scale = sqrt (major-axis * minor-axis)
lInner = 33.941125497
lOuter = 224.0
lAverage = 128.970562748
l = 57.6984848098
l = 81.4558441227
l = 105.213203436
l = 128.970562748
l = 152.727922061
l = 176.485281374
l = 200.242640687
ellipses B:
linear scale = sqrt (major-axis * minor-axis)
lInner = 32.0
lOuter = 112.0
lAverage = 72.0
l = 42.0
l = 52.0
l = 62.0
l = 72.0
l = 82.0
l = 92.0
l = 102.0
ellipses C:
linear scale = sqrt (major-axis * minor-axis)
lInner = 34.1320963318
lOuter = 167.567150719
lAverage = 100.849623525
l = 50.8114781301
l = 67.4908599285
l = 84.1702417269
l = 100.849623525
l = 117.529005324
l = 134.208387122
l = 150.88776892
ellipses D:
linear scale = sqrt (major-axis * minor-axis)
lInner = 40.4042077017
lOuter = 188.865031173
lAverage = 114.634619437
l = 58.9618106356
l = 77.5194135695
l = 96.0770165034
l = 114.634619437
l = 133.192222371
l = 151.749825305
l = 170.307428239

And here is the script itself:

# interpolate from enclosed ellipse to enclosing ellipse
# use 3x3 affine-transformation matrix to parameterize ellipse

import copy
import math

# ellipse-affine:  3x3 affine transformation matrix
# ellipse is given by affine transformation applied to the unit circle
# (x0, yo) is the center of the ellipse
# when the 2x2 linear transformation is a "pure rescaling" we have:
#
#   [ alpha,        beta / 2,  x0 ]
#   [ beta / 2,     gamma,     y0 ]
#   [ 0,            0,         1  ]
#
# if beta is zero, alpha and gamma are the semi-major and semi-minor axes


# some utility functions

# multiply two 3x3 list-of-list matrices together
def m3dmmm (m1, m2):
    retval = []
    for  i in range (3):
        retval.append ([0.0] * 3)
    for  i in range (3):
        for  j in range (3):
            for  k in range (3):
                retval[i][j] += m1[i][k] * m2[k][j]
    return  retval

# multiply a 3x3 list-of-list matrix onto a list vector
def m3dmvm (m, v):
    retval = [0.0] * 3
    for  i in range (3):
        for  j in range (3):
            retval[i] += m[i][j] * v[j]
    return  retval

# calculate the determinant of a 3x3 affine-transformation matrix
# that is, the determinant of its upper-left 2x2 linear-transformation sub-block
def a3ddet (a):
    return  a[0][0] * a[1][1]  -  a[1][0] * a[0][1]

# calculate the inverse of a 3x3 affine-transformation matrix
def a3dinv (a):
    d = a3ddet (a)
    retval = []
    for  i in range (3):
        retval.append ([0.0] * 3)
    retval[0][0]  =   a[1][1] / d
    retval[0][1]  =  -a[0][1] / d
    retval[1][0]  =  -a[1][0] / d
    retval[1][1]  =   a[0][0] / d
    retval[0][2]  =  -(retval[0][0] * a[0][2] + retval[0][1] * a[1][2])
    retval[1][2]  =  -(retval[1][0] * a[0][2] + retval[1][1] * a[1][2])
    retval[2][2]  =   1.0
    return  retval

# solve func (x, data) == 0 for x using bisection
# xlo, and xhi must bracket the root
def bisect (func, data, xlo, xhi):
    tol = 0.0
    ylo = func (xlo, data)
    if  ylo == 0.0:
        return  xlo
    yhi = func (xhi, data)
    if  yhi == 0.0:
        return  xhi
    if  (ylo > 0.0) == (yhi > 0.0):
        raise  ValueError ('xlo, xhi must bracket root')
    while (True):
        x = (xlo + xhi) / 2.0
        if  abs (x - xlo) <= tol  or  abs (x - xhi) <= tol:
            return  x
        y = func (x, data)
        if  y == 0.0:
            return  x
        if  (y > 0.0) == (ylo > 0.0):
            xlo = x
            ylo = y
        else:
            xhi = x
            yhi = y

# the symmetric part of the polar decomposition of a 2x2 matrix
# passed in as [m_00, m_01, m_10, m_11]
# symmetric 2x2 matrix is returned as [s_00, s_11, s_01] (where s_10 = s_01)
def polarsym2d (m):
    d = math.sqrt ((m[0] + m[3])**2 + (m[1] - m[2])**2)
    s0 = ( m[0] * (m[0] + m[3])  +  m[1] * (m[1] - m[2]) ) / d
    s1 = ( m[3] * (m[0] + m[3])  -  m[2] * (m[1] - m[2]) ) / d
    s2 = ( m[0] * m[2]  +  m[1] * m[3] ) / d
    return  [s0, s1, s2]

# "derotate" non-symmetric 2x2 linear block in a 3x3 affine transformation
# to produce equivalent ellipse with a symmetric ("pure rescaling") linear
# transformation.  that is, the derotated linear block is the symmetric
# part of the original linear block's polar decomposition.
def derotate (a):
    retval = copy.deepcopy (a)
    s = polarsym2d ([ a[0][0], a[0][1], a[1][0], a[1][1] ])
    retval[0][0] = s[0]
    retval[1][1] = s[1]
    retval[0][1] = s[2]
    retval[1][0] = s[2]
    return  retval


# returns "ellipse parameters" (from tuple) as affine transformation matrix
def eprmToEaff (ep):
    return  [
        [ ep[0],        ep[1] / 2.0,  ep[3] ],
        [ ep[1] / 2.0,  ep[2],        ep[4] ],
        [ 0.0,          0.0,          1.0   ]
    ]

# identity transformation is the unit circle
def getUnitCircle():
    return eprmToEaff ((1.0, 0.0, 1.0, 0.0, 0.0))

# transform the outer ellipse of a nested ellipse pair into the unit circle,
# and rotate so that the (transformed) inner ellipse is horizontal and is in
# the upper half plane
# returns the normalizing transformation and its inverse
def normalizePair (eo, ei):
    xfrm = a3dinv (eo)
    exi = m3dmmm (xfrm, ei)
    aa = exi[0][0]**2 + exi[1][0]**2
    bb = exi[0][0] * exi[0][1] + exi[1][0] * exi[1][1]
    dd = exi[0][1]**2 + exi[1][1]**2
    thtm = math.atan2 (2.0 * bb, aa - dd) / 2.0
    x = exi[0][0] * math.cos (thtm) + exi[0][1] * math.sin (thtm)
    y = exi[1][0] * math.cos (thtm) + exi[1][1] * math.sin (thtm)
    thte = math.atan2 (y, x)   # major-axis angle
    cx = exi[0][2]
    cy = exi[1][2]
    if  cx != 0.0  or  cy != 0.0:
        thtp = math.atan2 (exi[1][2], exi[0][2])   # center angle
        thtp -= thte   # rotated center angle
        # rotate center to upper half plane
        if  (thtp >= -math.pi  and  thtp < 0.0)  or  thtp >= math.pi:
            thte += math.pi
    rot = [
        [ math.cos (thte), math.sin (thte), 0.0],
        [-math.sin (thte), math.cos (thte), 0.0],
        [0.0, 0.0, 1.0]
    ]
    xfrm = m3dmmm (rot, xfrm)
    xfri = a3dinv (xfrm)
    return  (xfrm, xfri)

# transform ellipse from unit circle to guide circle
def  transform (e, eGuide):
    return  m3dmmm (eGuide, e)


# returns guide circle (as ellipse) for line segment
# (x, y) is the center of the line
# a is the half-length (semi-major axis)
def getGuideCircle (x, y, a):
    if  y < 0.0  or  y >= 1.0:
        raise  ValueError ('bad y')
    if  a <= 0.0  or  ((a + abs (x))**2 + y**2) > 1.0:
        raise  ValueError ('bad a')
    r0 = a / math.sqrt (1.0 - y**2)
    y0 = (1.0 - r0) * y
    return  [
        [  r0, 0.0,   x ],
        [ 0.0,  r0,  y0 ],
        [ 0.0, 0.0, 1.0 ]
    ]


# returns slope mismatch for ellipse and unit circle at height h
def getChordWrapper (h, a_b_y0):
    a, b, y0 = a_b_y0
    return  (1 - h**2) / h**2  -  b**2 * (b**2 - (h - y0)**2) / (a * (h - y0))**2

# finds the height of the "tangent chord" for horizontal ellipse in unit circle
def getChord (a, b, y0):
    eps = 1.e-5
    if  y0 <= 0.0:
        return 0.0
    if  y0 >= 1.0:
        return 1.0
    blo = (1.0 + eps) * y0
    bhi = (1.0 - eps) * y0 + b
    return  bisect (getChordWrapper, (a, b, y0), blo, bhi)

# returns ellipse tangent to circle at thetaT with "reflection angle" thetaE
def tangentEllipse (thetaE, thetaT):
    if  thetaE <= 0  or  thetaE >= thetaT:
        print  'thetaE =', thetaE, ', thetaT =', thetaT
        raise  ValueError ('bad thetaE')
    tht1 = thetaT + thetaE
    tht2 = thetaT - thetaE
    delta = 2.0 * math.cos (thetaT) / (1.0 / math.tan (tht1) + 1.0 / math.tan (tht2))
    a = (delta / 2.0) * (1.0 / math.sin (tht1) + 1.0 / math.sin (tht2))
    d1 = math.cos (thetaT) - delta / math.tan (tht1)
    d2 = delta / math.tan (tht2) - math.cos (thetaT)
    rdiff = 2.0 * (d2 - d1) / (d1 + d2)
    if  abs (rdiff) > 1.e-6:
        print  'd1 =', d1, ', d2 =', d2, ', rdiff =', rdiff
        raise  RuntimeError ('d1, d2 mismatch')
    c = d1
    b = math.sqrt (a**2 - c**2)
    y = math.sin (thetaT) - delta
    e = [
        [   a, 0.0, 0.0 ],
        [ 0.0,   b,   y ],
        [ 0.0, 0.0, 1.0 ]
    ]
    return e

# returns  b (thetaE, thetaT) - bTarg
def tangentEllipseBWrapper (thetaE, thetaT_bTarg):
    thetaT, bTarg = thetaT_bTarg
    e = tangentEllipse (thetaE, thetaT)
    return  e[1][1] - bTarg

# returns tangent ellipse with semi-minor axis b
def tangentEllipseByB (b, thetaT):
    eps = 1.e-5
    blo = eps * thetaT
    bhi = (1.0 - eps) * thetaT
    tol = (eps**2) * thetaT
    thetaE = bisect (tangentEllipseBWrapper, (thetaT, b), blo, bhi)
    return  tangentEllipse (thetaE, thetaT)


# returns  area (thetaE, thetaT) - aTarg
def tangentEllipseAreaWrapper (thetaE, thetaT_aTarg):
    thetaT, aTarg = thetaT_aTarg
    e = tangentEllipse (thetaE, thetaT)
    return  e[0][0] * e[1][1] - aTarg

# returns tangent ellipse with "area" (semi-major * semi-minor) a
def tangentEllipseByArea (a, thetaT):
    eps = 1.e-5
    blo = eps * thetaT
    bhi = (1.0 - eps) * thetaT
    thetaE = bisect (tangentEllipseAreaWrapper, (thetaT, a), blo, bhi)
    return  tangentEllipse (thetaE, thetaT)


# returns (a, b) proportionally scaled between (a0, an) and (b0, bn)
# such that  a * b = r
def linearAreaStepFactors (a0, b0, an, bn, r):
    r0 = a0 * b0
    rn = an * bn
    abx = a0 * bn + an * b0
    qa = r0 + rn - abx
    qb = abx - 2.0 * r0
    qc = r0 - r
    eta = ( -qb + math.sqrt (qb**2 - 4.0 * qa * qc) ) / (2.0 * qa)
    return  ( (1.0 - eta) * a0 + eta * an, (1.0 - eta) * b0 + eta * bn )

# interpolate ellipse to unit circle with n intermediate steps of
# equally-spaced "linear scale" (square-root-of-area)
# starting ellipse is given by "tangent-ellipse" line segment and semi-minor axis b0
# returns lists of interpolating ellipses, tangent ellipses, and guide circles
def interpolateEllipse (x, y, a, n, b0):
    guideCircle = getGuideCircle (x, y, a)
    r0 = guideCircle[0][0]
    thetaT = math.asin (y)
    br = b0 / r0   # initial semi-minor axis, b0, rescaled to unit circle
    # get area of initial (b0) ellipse
    a0 = 0.0
    if  br != 0.0:
        a0 = a3ddet (tangentEllipseByB (br, thetaT))
    # rescale root-area to guide-circle radius
    a0r = r0 * math.sqrt (a0)
    #
    ies = []
    tes = []
    gcs = [guideCircle]
    for  i in range (0, n):
        step = (i + 1.0) / (n + 1.0)
        # ar is the linear scale (root-area) desired for each interpolating ellipse
        ar = a0r + step * (1.0 - a0r)   # ar ranges from a0r to 1
        art, gCR = linearAreaStepFactors (a0r / r0, r0, 1.0, 1.0, ar)
        at = art**2  # desired area of tangent ellipse
        te = tangentEllipseByArea (at, thetaT)
        # interpolate guide-circle center (x and y) down to 0
        gCX = ((1.0 - gCR) / (1.0 - r0)) * guideCircle[0][2]
        gCY = ((1.0 - gCR) / (1.0 - r0)) * guideCircle[1][2]
        gc = eprmToEaff ((gCR, 0, gCR, gCX, gCY))
        # transform tangent ellipse to guided interpolating ellipse
        ie = transform (te, gc)
        ies.append (ie)
        tes.append (te)
        gcs.append (gc)
    gcs.append (getUnitCircle())
    return (ies, tes, gcs)


# interpolate from inner ellipse to outer ellipse with n intermediate steps
def interpolateEllipses (ei, eo, n):
    xfrm, xfri = normalizePair (eo, ei)
    # eix is horizontal ellipse within unit circle with center in upper half plane
    eix = derotate (m3dmmm (xfrm, ei))
    a = eix[0][0]   # semi-major axis
    b = eix[1][1]   # semi-minor axis
    x = eix[0][2]   # center (x)
    y = eix[1][2]   # center (y)
    yc = getChord (a, b, y)  # center (y) of tangent chord
    ac = a * math.sqrt (1.0 - ((yc - y) / b)**2.0)
    iex = interpolateEllipse (x, yc, ac, n, b)[0]
    iexi = []
    for  e in iex:
        iexi.append (derotate (m3dmmm (xfri, e)))
    return iexi


# imagej-specific code

from ij import IJ
from ij.gui import EllipseRoi
from ij.gui import Overlay
from java.awt import Color

def createCanvas (title = 'ellipse'):
    imp = IJ.createImage (title, '8-bit ramp', 512, 512, 1)
    imp.getProcessor().multiply (0.25)
    imp.getProcessor().add (32.0)
    imp.setOverlay (Overlay())
    return  imp

def drawRoi (imp, roi, width, color):
    roi.setStrokeWidth (width)
    roi.setStrokeColor (color)
    imp.getOverlay().add (roi)

def eroiToEaff (eroi):
    params = eroi.getParams()
    x1 = params[0]
    y1 = params[1]
    x2 = params[2]
    y2 = params[3]
    ar = params[4]
    # major axis, minor axis, major-axis angle
    x0 = (x1 + x2) / 2.0
    y0 = (y1 + y2) / 2.0
    dx = x2 - x1
    dy = y2 - y1
    a = (dx**2 + dy**2)**0.5 / 2.0
    b = ar * a
    # if  dx <= 0.0:
    #     tht = math.pi / 2.0
    # else:
    #     tht = math.atan (dy / dx)
    tht = math.atan2 (dy, dx)
    # ellipse-affine parameters
    c = math.cos (tht)
    s = math.sin (tht)
    alp = c**2 * a + s**2 * b
    bth = s * c * (a - b)
    gam = s**2 * a + c**2 * b
    #
    return  [
        [ alp, bth,  x0 ],
        [ bth, gam,  y0 ],
        [ 0.0, 0.0, 1.0 ]
    ]

def eaffToEroi (eaff):
    alp = eaff[0][0]
    bth = eaff[0][1]
    gam = eaff[1][1]
    x0 = eaff[0][2]
    y0 = eaff[1][2]
    # verify "symmetric" ellipse-affine
    tol1 = 1.e-12
    tol2 = 1.e-6
    bth2 = eaff[1][0]
    if  (abs (bth) > tol1  or  abs (bth) > tol1)  and  (abs (bth2 - bth) / (abs (bth2) + abs (bth))) > tol2:
        raise ValueError ('eaffToEroi: eaff not symmetric')
    det = alp * gam - bth**2
    if  det <= 0:
        raise ValueError ('eaffToEroi: non-positive det')
    if  eaff[2][0] != 0.0  or  eaff[2][1] != 0.0  or  eaff[2][2] != 1.0:
        raise ValueError ('eaffToEroi: eaff not affine')
    # major axis, minor axis, major-axis angle
    tr = alp + gam
    dab = math.sqrt (tr**2 - 4.0 * det)
    a = (tr + dab) / 2.0
    b = (tr - dab) / 2.0
    tht = math.atan2 (2.0 * bth, alp - gam) / 2.0
    # EllipseRoi
    x1 = x0 + math.cos (tht) * a
    y1 = y0 + math.sin (tht) * a
    x2 = x0 - math.cos (tht) * a
    y2 = y0 - math.sin (tht) * a
    ar = b / a
    eroi = EllipseRoi (x1, y1, x2, y2, ar)
    return  eroi


# does one Roi contain another
def  roiContains (rOuter, rInner):
    for  p in rInner:
	if  not rOuter.contains (p.x, p.y):
	    return False
    return  True

# version
print  IJ.getFullVersion()

# example interpolations
eis = []
eos = []
eis.append (EllipseRoi (255.0,  63.0, 255.0, 159.0, 0.50))
eos.append (EllipseRoi  (31.0, 255.0, 479.0, 255.0, 1.00))
eis.append (EllipseRoi  (95.0, 255.0, 159.0, 255.0, 1.00))
eos.append (EllipseRoi  (31.0, 255.0, 479.0, 255.0, 0.25))
eis.append (EllipseRoi (300.0, 400.0, 325.0, 295.0, 0.40))
eos.append (EllipseRoi (345.0, 450.0, 180.0,  85.0, 0.70))
eis.append (EllipseRoi  (40.0, 305.0, 150.0, 370.0, 0.40))
eos.append (EllipseRoi  (40.0, 380.0, 470.0, 150.0, 0.60))

titles = ['ellipses A', 'ellipses B', 'ellipses C', 'ellipses D']
thick = 3
widths = [0, 0, 0, thick, 0, 0, 0]
colors = [
    Color.ORANGE,
    Color.YELLOW,
    Color.GREEN,
    Color.BLUE,
    Color.GREEN,
    Color.YELLOW,
    Color.ORANGE
]
for  ttl, ei, eo in zip (titles, eis, eos):
    print  ttl + ':'
    if  not roiContains (eo, ei):
        print  'ei not contained within eo, skipping'
        continue
    # check linear scale (root area) of ellipses
    print  'linear scale = sqrt (major-axis * minor-axis)'
    lInner = math.sqrt (a3ddet (eroiToEaff (ei)))
    lOuter = math.sqrt (a3ddet (eroiToEaff (eo)))
    lAverage = (lInner + lOuter) / 2.0
    print  'lInner =', lInner
    print  'lOuter =', lOuter
    print  'lAverage =', lAverage
    img = createCanvas (ttl)
    drawRoi (img, ei, thick, Color.RED)
    drawRoi (img, eo, thick, Color.RED)
    eiaff = eroiToEaff (ei)
    eoaff = eroiToEaff (eo)
    ies = interpolateEllipses (eiaff, eoaff, 7)
    for  e, w, c in zip (ies, widths, colors):
        print  'l =', math.sqrt (a3ddet (e))
        eroi = eaffToEroi (e)
        drawRoi (img, eroi, w, c)
    # save as png's for upload
    imgFlat = img.flatten()
    IJ.saveAs (imgFlat, 'PNG', ttl)
    imgFlat.close()
    img.show()

Thanks, mm

5 Likes