Fitting a circle to maximum intensity in 2D data

Question

I have 2D intensity data that looks a bit like a volcano with a crater. I would like to fit an ellipse to the rim of the volcano defined by (smoothed) maximum intensity at each angle around the center of the image. Is this a problem for which an algorithm exists in python or openCV or should I write my own?

I was thinking to smooth the image and then find the maximum along a number of radial profiles and then least-squares-ellipse fitting such as this: https://github.com/bdhammel/least-squares-ellipse-fitting

but maybe something else already exists?

There's [circle and ellipse Hough Transforms](https://scikit-image.org/docs/dev/auto_examples/edges/plot_circular_elliptical_hough_transform.html) but they mostly look for one pixel wide. Might want to start there though. — Daniel F, Jul 19 '21 at 10:03
I would try to fit a donut-like distribution function like a difference of two weighted gaussinas with the same mean and different variances. — yann ziselman, Jul 19 '21 at 11:21
Try to determine the center. Then do a cartesian to polar transform. Then fit a line to the bright part since the ring will be along a line in the polar transform. See cv2.linearPolar() at https://docs.opencv.org/4.1.1/da/d54/group__imgproc__transform.html#gaa38a6884ac8b6e0b9bed47939b5362f3 and example at https://stackoverflow.com/questions/51675940/converting-an-image-from-cartesian-to-polar-limb-darkening — fmw42, Jul 19 '21 at 15:33

score 0 · Answer 1 · answered Jul 19 '21 at 23:12

I am not quite sure what your desired output is. For now I'll assume you want to know the center and radius of your volcano.

My approach is to set up a volcano generator and wiggle on it's parameters till it generates basically yours. Then since I know how it was generated I can tell you all about it.

import torch
import torch.nn as nn
import torch.optim

class Net(nn.Module):
    def __init__(self):
        super().__init__()
        self.foci = nn.Parameter(torch.tensor([[30.0,30.0],[30.0,30.0]]))
        self.r = nn.Parameter(torch.tensor((20.0)))
        self.a = nn.Parameter(torch.tensor((-0.1)))
        self.b = nn.Parameter(torch.tensor((0.01)))

    def generate_volcano(self):
        eps = 1e-15
        x,y = torch.meshgrid(torch.arange(60.0), torch.arange(60.0))
        d1 = ((self.foci[0,0] - x)**2+(self.foci[0,1] - y)**2 + eps)**(1/2)
        d2 = ((self.foci[1,0] - x)**2+(self.foci[1,1] - y)**2 + eps)**(1/2)
        return torch.sigmoid(self.a*(d1+d2-self.r)**2+self.b)

n = 10**5
optimizer = optim.Adam(net.parameters())
torch_image = torch.from_numpy(image.astype('float32'))
loss_funcition = nn.MSELoss()

for i in range(n):
    optimizer.zero_grad()
    volcano = net.generate_volcano()
    loss = loss_funcition(volcano, torch_image)
    loss.backward()
    optimizer.step()
    if i % (n//8) == 0:
        print('loss', loss.item())

This gives me a "volcano" like this

with focal points [32.1, 29.6],[32.1, 29.6] and radius 38.527.

Fitting a circle to maximum intensity in 2D data

1 Answers1