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