Skip to content

Instantly share code, notes, and snippets.

@teoliphant
Last active May 16, 2022 01:30
Show Gist options
  • Save teoliphant/7709098 to your computer and use it in GitHub Desktop.
Save teoliphant/7709098 to your computer and use it in GitHub Desktop.
Create masks (circle and ellipse) starting with a simple level cutoff.
import numpy as np
from scipy.optimize import fmin
from math import sqrt
import sys
def circ(h, k, r, x, y):
return (x-h)**2 + (y-k)**2 <= r**2
def func((h,k,r), x,y, mask):
return (circ(h,k,r,x,y) - mask).sum()
def moments(mask, x, y):
N = float(mask.sum())
h0 = (x*mask).sum() / N
k0 = (y*mask).sum() / N
r0 = sqrt(2*((((x-h0)**2 + (y-k0)**2)*mask).sum() / N))
return h0, k0, r0
def update_mask(mask, mask0, x, y):
h0, k0, r0 = moments(mask, x, y)
return mask0 * circ(h0, k0, r0, x, y)
def circle_mask(img, level=0.65, iter=8):
mask0 = img > level * img.mean()
x, y = np.indices(img.shape)
mask = mask0
for k in range(iter):
mask = update_mask(mask, mask0, x, y)
h0, k0, r0 = moments(mask, x, y)
print "Starting values...", (h0, k0, r0)
hopt, kopt, ropt = fmin(func, (h0, k0, r0), args=(x, y, mask))
print "Ending values...", (hopt, kopt, ropt)
return circ(hopt, kopt, ropt, x, y)
def main():
import scipy.misc as sm
image = sm.imread(sys.argv[1], flatten=True)
circmask = circle_mask(image)
sm.imsave('circ_mask.png', circmask)
sm.imsave('circ_masked.png', circmask*image)
if __name__ == '__main__':
main()
import numpy as np
from scipy.optimize import fmin
from math import sqrt, asin, sin, cos
import sys
def ellipse(h, k, a, b, phi, x, y):
xp = (x-h)*cos(phi) + (y-k)*sin(phi)
yp = -(x-h)*sin(phi) + (y-k)*cos(phi)
return (xp/a)**2 + (yp/b)**2 <= 1
def func((h,k,a,b,phi), x,y, mask):
return (ellipse(h,k,a,b,phi,x,y) - mask).sum()
def sgn(x):
return -1.0*(x<0) + 1.0*(x>0)
def moments(mask, x, y):
N = float(mask.sum())
h0 = (x*mask).sum() / N
k0 = (y*mask).sum() / N
sxx = (((x-h0)**2)*mask).sum() / N
syy = (((y-k0)**2)*mask).sum() / N
sxy = (((x-h0)*(y-k0))*mask).sum() / N
term1 = 2*(sxx + syy)
term2 = sgn(sxy) * sqrt(4*(sxx-syy)**2 + 16*sxy*sxy)
asq = term1 + term2
bsq = term1 - term2
diff = asq - bsq
if diff == 0:
phi = 0.0
else:
phi = 0.5*asin(8*sxy / diff)
return h0, k0, sqrt(asq), sqrt(bsq), phi
def update_mask(mask, mask0, x, y):
h0, k0, a0, b0, phi0 = moments(mask, x, y)
return mask0 * ellipse(h0, k0, a0, b0, phi0, x, y)
def ellipse_mask(img, level=0.65, iter=8):
mask0 = img > level * img.mean()
x, y = np.indices(img.shape)
mask = mask0
for k in range(iter):
mask = update_mask(mask, mask0, x, y)
h0, k0, a0, b0, phi0 = moments(mask, x, y)
print "Starting values...", (h0, k0, a0, b0, phi0)
hopt, kopt, aopt, bopt, phiopt = fmin(func, (h0, k0, a0, b0, phi0), args=(x, y, mask))
print "Ending values...", (hopt, kopt, aopt, bopt, phiopt)
return ellipse(hopt, kopt, aopt, bopt, phiopt, x, y)
def main():
import scipy.misc as sm
image = sm.imread(sys.argv[1], flatten=True)
ellmask = ellipse_mask(image)
sm.imsave('ell_mask.png', ellmask)
sm.imsave('ell_masked.png', ellmask*image)
if __name__ == '__main__':
main()
@arpit1ug
Copy link

Hey,

Can you please give me some reference or a specific name to learn the algorithm used in moments

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment