[Python] subtract parabolic average from image (vignette removal)

I have to process the images I acquire with the microscope with a vertical binning, to obtain the distances between the dark/bright regions. Images are like this:

I defined a threshold function that crops the circular image, makes the binning, and subtract a parabolic profile:

def threshold(image, X0Y0R):
'''
crop the image along the circle of the ocular (step1)
vertical binning of image, and subtract parabolic background (step2),
resulting in sub_binn

X0Y0R is the (X0,Y0,R) tuple, with circle center and radius

returns:
((np.max(sub_binn)-np.min(sub_binn))/2)+np.min(sub_binn) -> sub_binn's avg gray level
image -> the rotated, circle-cropped image
binn -> 1D array vert binning, previous to parabolic subtraction
sub_binn -> 1D array used by image_analisys()
'''
image = rotate(image, tilt, reshape=True)

######## -->> step1
X0,Y0,R = X0Y0R[0],X0Y0R[1],X0Y0R[2]
Y, X = np.ogrid[:image.shape[0], :image.shape[1]]
outer_disk_mask = (X - X0)**2 + (Y - Y0)**2 > R**2
image[outer_disk_mask] = False

######## -->> step2
binn = []
for i in range(len(image[0])):
binn = np.append(binn, np.sum(np.transpose(image)))
binn = binn[np.argwhere(binn>0)]
x_n = np.arange(len(binn))
par = np.polyfit(range(len(binn)), binn, 2)

sub_binn = []
for i in range(len(binn)):
sub_binn = np.append(sub_binn, binn[]-(x_n*par[1]+par[0]*(x_n*x_n)))

return ((np.max(sub_binn)-np.min(sub_binn))/2)+np.min(sub_binn),image, binn, sub_binn


Now when I examine the result:

image = img.imread('image.png').astype(np.int32)
print(image.shape)
X0,Y0,R = 950,700,550
it=threshold(image,X0Y0R=(X0,Y0,R))
fig,ax=plt.subplots(1,2,figsize=(12,4))
plt.subplot(121)
plt.imshow(it[1], interpolation='none',cmap="gray")
plt.subplot(122)
plt.plot(it[2],color='gray',label='vertical binning')
plt.plot(it[3],color='black',label='after parab subtract')
plt.hlines(y=it[0], color='red',xmin=0,xmax=1000)
print(it[1].shape, 2*R)
plt.legend()
plt.show()


I see that the parabolic fit is never correct - e.g. in this case it overestimates it:



Since it will have to work on different images, there is some way to extract a reasonably horizontal profile, unsupervised?

Continue reading...