Finding an average representation for a set of 2D contours

Question

I work at a lab and use mice for behavioral studies. Part of what I want to show is how a known stimulus will elicit movement of the ear in a predictable way. I have a script to segment an ear across e.x. 21 video frames and generate contours.

When the contours are overlayed and coded [RED = before stimulus; BLUE = after stimulus], one can appreciate an overall difference in the ear shape/orientation.

I'd like to display this visual data in a more friendly way. My idea is to find an 'average' representation for the RED and BLUE contours and compare the two instead of this messy business with 21 contours. Bonus would be if, from this 'average' contour, we could also display a standard deviation of the inherent variance within the RED or BLUE contour sets.

To elucidate further, I want to relate this challenge to the task of visualizing 1D time series data. For example, figure 4A of this paper visualizes time series data as shown below. Many neural fluorescence recordings are averaged to a single green line and the variance is displayed as a gray background. I'm wondering if something similar is possible in 2D space with my contours?

The main challenge is that, compared to time series data, there isn't a 1-to-1 correspondence between a coordinate from one contour to a coordinate of another contour. For 1D data, we simply average all values for some t to find the mean & variance at that spot, then apply the operation over the entire domain. There is no index t that correlates coordinate pairs between contours. Furthermore, contours won't even have the same number of coordinates due to their difference in size/shape.

But what if there is a rough index t for coordinates in a 2D contour. Suppose that we calculate the center of mass (CoM) for the BLUE contours? From the CoM, perhaps we could iterate over small sectors of a very large circle (r=1000, d_theta=1) that encapsulates the contours. This sector acts like a reference t, so all coordinates within this sector will be averaged to a single coordinate. If we iterate through all sectors of this circle, then I assume we'd arrive to a fairly good representation of a 'mean' 2D contour. Additionally, we could obtain a measure of variance from each sector and plot it as a gray bar that extends tangentially from each new 'average' point of the new 'average' contour. The issue is that those bars likely won't look very nice and smooth like in the 1D time series case. Below is a rough visualization of the output I'm envisioning.

Please let me know if any ideas sparked in your head! I leave sample code below to generate roughly aligned circular contours if anyone wants to play around/test things.

import numpy as np
import matplotlib.pyplot as plt

#generate some ovals
num_ovals = 4
points = []

disp = np.random.rand(num_ovals, 2) * 0.5 # random displacements

for _ in range(num_ovals):
    center = np.random.rand(2) * 0.05
    d_x = 1 + np.random.rand() * 0.5  # Dilation noise x
    d_y = 1 + np.random.rand() * 0.5  # Dilation noise y 
    angle = np.random.rand() * 0.2  # Angle offset
    
    # generate with random dilation
    t = np.linspace(0, 2 * np.pi, 100)
    x = center[0] + d_x * 2 * np.cos(t)
    y = center[1] + d_y * np.sin(t)
    
    # apply random displacement
    x += disp[_][0]
    y += disp[_][1]
    
    # apply random rotations
    rotation_matrix = np.array([[np.cos(angle), -np.sin(angle)],
                                [np.sin(angle), np.cos(angle)]])
    points.append(np.dot(np.vstack((x, y)).T - center, rotation_matrix) + center)
    
comX, comY = np.mean(np.vstack(points), axis=0).T #center of mass

#plot figure
plt.figure(figsize=(8, 8))

for p in points:
    plt.plot(p[:, 0], p[:, 1], marker='.', markersize = 2)

plt.plot(comX, comY, marker='x', markersize=10, color='red', label='Center of Mass')
    
plt.legend()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
 

Perhaps finding an average contour in the same sense that time series data is averaged may not be feasible. In that case, I propose quickly another thought/curiosity of mine. Lets forget about the lines connecting all the coordinate pairs that make the contours, and just think of the contours as a collection of coordinate pairs that form the 2D point cloud. Would it be possible to find single contour that best approximates this 2d point cloud? I think this would be a lot easier to implement even thought I don't see any way to reasonably infer about variance in a simple way.

Thanks for reading if you've made it through!

Answer 1

After hearing about the concept of Wasserstein barycenters introduced in the comments, and being introduced to the leading python library for optimal transport algorithms (POT), I have found the current optimal solution to my question.

Here is an example script from the POT library that describes how to calculate a weighted Sinkhorn barycenter for two discrete distributions. A sample output from the linked script is shown below. The corners represent discrete 2D distributions and the sides display the weighted interpolated barycenters.

I've adapted this script to work for my own purposes. One thing to keep in mind with the result displayed above is that all 2D distributions have the same number of samples n_samples = 200. Additionally, you can see from the plots that all input distributions are centered at 0 and (x,y) ∈ [-3, 3]. Due to my unfamiliarity with these optimal transport algorithms, I was not initially aware that these features of the distributions were important for the algorithms to work correctly.

In the area of image segmentation, masks and their contours are typically positioned far away from the origin (0,0) corner of the image. Additionally, the contour outputs from, for example, openCV's findContours() function will have various shapes and sizes. Here is my own segmentation of a mouse's ear across 20 frames:

The coordinate arrays representing the contours above range in length from 486 to 572 elements. Hence, I first wrote the following scripts to resample the contours for equal coordinate lengths and then transform them to the fit in a 3x3 square centered at the origin.

from scipy.interpolate import interp1d
import numpy as np

def sampler(cnts, sample_size):
    
    interp_cnts = []
    for c in cnts:
        x,y = c.T

        interp_x = interp1d(np.linspace(0, 1, len(x)), x, kind='cubic')
        interp_y = interp1d(np.linspace(0, 1, len(y)), y, kind='cubic')

        # Interpolate to create new x and y arrays with 'maxN' points
        new_x = interp_x(np.linspace(0, 1, sample_size))
        new_y = interp_y(np.linspace(0, 1, sample_size))

        # Combine the new x and y arrays to form the interpolated contour
        interp_cnt = np.column_stack((new_x, new_y))
        interp_cnts.append(interp_cnt)
        
    return interp_cnts

n_samples = 200
s_cnts = sampler(cnts, n_samples)

After resampling all contours to hold 200 elements (which did not affect the accuracy of the contours, I transformed them to fit the origin).

def c_transform(cnts, d): #list of contours, and dilation factor post-normalization (typically 10)
    
    combined_points = np.vstack(cnts)
    common_centroid = np.mean(combined_points, axis=0)

    # normalization factor
    max_abs = max(np.max(np.abs(combined_points[:, 0])), np.max(np.abs(combined_points[:, 1])))
    
    transformed_cnts = []
    for c in cnts:
        c_n = (c - common_centroid) / max_abs * d  #center mass & normalize
        transformed_cnts.append(c_n)

    return transformed_cnts

t_cnts = c_transform(s_cnts, 10)

Below you can view the output:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(3,3))
ax.invert_yaxis()

for c in t_cnts:
    ax.scatter(x=c[:, 0], y=c[:, 1], c='steelblue', s=1)

Next, I finally implemented the script to calculate the Sinkhorn barycenter, which was basically a copy/paste from the link above.

import ot

def Sinkhorn_OT(cnts, n_samples):
    
    reg = 1e-2  # Entropic Regularization
    numItermax = 20  # Maximum number of iterations for the Barycenter algorithm (20)
    numInnerItermax = 50  # Maximum number of sinkhorn iterations (50)
    
    N = len(cnts)
    weights = np.full(N, 1 / N).astype(np.float32) #give equal weight to all cnts
    
    a = ot.unif(len(cnts[0]))
    measures_weights=[a] * N  #provide uniform histogram for weights
    

    
    XB_init = np.random.randn(n_samples, 2)
    XB = ot.bregman.free_support_sinkhorn_barycenter(
        measures_locations=cnts,
        measures_weights=measures_weights,
        weights=weights,
        X_init=XB_init,
        reg=reg,
        numItermax=numItermax,
        numInnerItermax=numInnerItermax
    )
    
    return XB

This function takes a list of transformed contours and outputs the barycenter. Specifically, all the contours are given equal weight so that the output is an 'average' of all contours in the list.

test_cnts = [t_cnts[0], t_cnts[19]]
XB = Sinkhorn_OT(test_cnts, n_samples)

fig, ax = plt.subplots(figsize=(3,3))
ax.invert_yaxis()

for c in test_cnts:
    ax.plot(c[:, 0], c[:, 1], c='steelblue', alpha=0.5)
    
ax.scatter(x=XB[:, 0], y=XB[:, 1], c='gold', s=3)

Below is the output, with the input distributions shown as blue lineplots and the barycenter outlined as a yellow scatterplot.

I wrote the code above to handle more than just two distributions, because I want to find the average contour across an arbitrary number of frames. Here is another output with 10 frames for input:

I believe these results fit any modest criteria of an 'average' shape, which is why I will accept this as the solution. Additionally, one may want to invert the initial transformations so the barycenter fits on the original image from which the input contours were extracted. That is a simple computation which can be done by reversing the c_transform() function I defined above.

For metrics, it takes me ~ 4s to calculate the Sinkhorn barycenter for 10 distributions of 200 samples each. I'm using an NVIDIA RTX A2000 12GB with pytorch. I have not parallelized anything, so there is immense room for improvement if people are working with larger datasets.

I apologize if my code is crude or lacks convention -- I do not have any formal education/background in computer programming. And I want to thank user Corralien for their advice once again.

Edit: One issue I'm having, which is very minor, is that the resulting barycenter can only be plotted as a scatterplot. It appears that the coordinates that make the average' outline are disordered, so a line-plot looks super jumbled. I don't know if there's a simple way to reorder the points. The only thing that could work is a convexhull on the set of 2d coordinates that form the barycenter, but the convexhull will lose a lot of detail from the contour.

Edit II: I have developed logic to organize the disorganized contour points such that if you plot them using plt.plot() they form a clean contour again. I won't post it here because it has more to do with aesthetics rather than the topic of the question, but if you'd like to see it let me know.