Create masks (circle and ellipse) starting with a simple level cutoff.

circ_mask.png

circ_masked.png

circle_mask.py

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()

ell_mask.png

ell_masked.png

ellipse_mask.py

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()

myimage.png

Hey,

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