Calculate Image Moments of an Ellipse

Image moments can be used to determine the center of gravity or average pixel location of a binary object as well as the orientation of an object according to its axes of variation.

The interactive plot below demonstrates how to calculate the image moments of an ellipse-defined mask.

import itk
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import interact

output_size = [100, 100] # size of binary image containing the ellipse

@interact(tx=(0,output_size[0],1), ty=(0,output_size[1],1), a=(5,50,1), b=(5,50,1), theta=(0,2*np.pi,0.1))
def foo(tx=30, ty=50, a=5, b=10, theta=np.pi/4.0):
    '''
    Creates a binary image of an ellipse and calculates the image moments.  Major (purple) and minor (green)
    principal axes are displayed.

    Parameters
    ==========
    tx, ty : translation x and y
    a, b : ellipse horizontal and vertical widths before rotation
    theta : angle of rotation (radians)
    '''

    # ellipse starts as unit circle, use world transform to define final ellipse
    ellipse = itk.EllipseSpatialObject[2].New()
    ellipse_transform = itk.AffineTransform[itk.D, 2].New()
    ellipse_transform.Scale([a, b])
    ellipse_transform.Rotate2D(theta)
    ellipse_transform.Translate([tx, ty])
    ellipse.SetObjectToWorldTransform(ellipse_transform)

    ellipse_img = itk.spatial_object_to_image_filter(input=ellipse, inside_value=1, outside_value=0, size=output_size)

    momentor = itk.ImageMomentsCalculator.New(Image=ellipse_img)
    momentor.Compute()

    centroid = momentor.GetCenterOfGravity()
    prin_axes = itk.array_from_matrix(momentor.GetPrincipalAxes())
    minor_axes = prin_axes[0]
    major_axes = prin_axes[1]

    fig, ax = plt.subplots(figsize=[8,8])
    plt.imshow(ellipse_img, cmap='gray')
    plt.scatter(centroid[0], centroid[1])

    minor_pt = centroid + minor_axes*5
    plt.plot([centroid[0], minor_pt[0]], [centroid[1], minor_pt[1]], color='green')

    major_pt = centroid + major_axes*5
    plt.plot([centroid[0], major_pt[0]], [centroid[1], major_pt[1]], color='purple')

    print(momentor)