I have a image converted to a numpy array. The image is of 2 curved lines. I'm trying to get the coordinates of the point where the highest black line peaks. I want to put a red square around this to highlight the point. So far my code only includes the image import and I'm struggling to find the right way to do this. Is it done through using slice, transpose, max and min functions? I tried these by they end up not getting me any closer. The line is quite thick so it should only select at the top most value of the picture.
import numpy as np
import os
import cv2 as cv2
b = cv2.imread("graph.png",0)
cv2.imshow("frame",b)
cv2.waitKey(0)
cv2.destroyAllWindows()
CodePudding user response：
I made a solution that uses linear least squares fitting an elipse to your image. (the two curves in your example clearly look like an elipse). This could of course be modified with the same approach for curves of different shapes.
We will need the following packages:
import os
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.optimize
We have to load in the png image to memory:
ImageData = plt.imread(os.path.expanduser("~") "/Downloads/ImageCurveFindMax.png" )
plt.figure()
plt.imshow(ImageData)
plt.draw()
Turn the image data into coordinates
white_rgba = np.ones(4)
black_rgba = np.array([0,0,0,1])
LineCoords = []
for y, ImageData_row in enumerate(ImageData):
for x, ImageData_pixel in enumerate(ImageData_row):
if not np.array_equal( ImageData_pixel, white_rgba ) :
#print ( ImageData_pixel )
LineCoords.append([x,y])
LineCoords = np.array(LineCoords)
X = np.atleast_2d( LineCoords[:,0] ).T
Y = np.atleast_2d( LineCoords[:,1] ).T
xmin = np.min(X)
xmax = np.max(X)
xrng = xmax - xmin
ymin = np.min(Y)
ymax = np.max(Y)
yrng = ymax - ymin
Get the maximum y-coords naively (this will be off-center):
ymaxnaiveindex = np.argmax( LineCoords[:, 1] )
ymaxnaivecoord = LineCoords[ymaxnaiveindex]
Fit the curves to fully general elipse using linear least squares: Taken directly from: How to fit a 2D ellipse to given points
# Ax^2 Bxy Cy^2 Dx Ey = 1
# Formulate and solve the least squares problem ||Ax - b ||^2
A = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(A, b)[0].squeeze()
# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 {1:.3}xy {2:.3}y^2 {3:.3}x {4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))
Get the coordinate on the fit ellipse such that y is maximum. Solve for y(x) and maximize:
param_A, param_B, param_C, param_D, param_E = x
# Cy^2 (E Bx)y (Ax^2 Dx - 1) = 0
# y = [ -(E Bx) - sqrt( (E Bx)^2 - 4 C (Ax^2 Dx - 1) ) ] / 2C
# quadratic_numerator = -(param_E param_B*x) - np.sqrt( (param_E param_B*x)**2 - 4 * param_C * (param_A*x**2 param_D*x - 1) )
# quadratic_denominator = 2*param_C
def yplus(x):
quadratic_a = param_C
quadratic_b = param_E param_B*x
quadratic_c = param_A*x**2 param_D*x - 1
quadratic_numerator = -quadratic_b np.sqrt( quadratic_b**2 - 4*quadratic_a*quadratic_c)
quadratic_denominator = 2*quadratic_a
return quadratic_numerator / quadratic_denominator
def yminus(x):
quadratic_a = param_C
quadratic_b = param_E param_B*x
quadratic_c = param_A*x**2 param_D*x - 1
quadratic_numerator = -quadratic_b - np.sqrt( quadratic_b**2 - 4*quadratic_a*quadratic_c)
quadratic_denominator = 2*quadratic_a
return quadratic_numerator / quadratic_denominator
#We could use algebra to find the max / min here - but we can be lazier and use: scipy.optimize.minimize
# y = yplus(x) OR yminus(x)
# ymax @ d/dx yplus(x) == 0 <--- solve this eqn is tedious for this case
yplusmaxObject = scipy.optimize.minimize(
fun = lambda a: -yplus(a),
x0 = xmin xrng/2,
bounds = (xmin, xmax ),
method = 'Nelder-Mead',
)
yplusmaxY = -yplusmaxObject.fun
yplusmaxX = yplusmaxObject.x
#It turns out after algebra the y = -b - sqrt(...) soln will be the upper elipse. We don't need the yplus code....
yminusmaxObject = scipy.optimize.minimize(
fun = lambda a: -yminus(a),
x0 = xmin xrng/2,
bounds = (xmin, xmax ),
method = 'Nelder-Mead',
)
yminusmaxY = -yminusmaxObject.fun
yminusmaxX = yminusmaxObject.x
Construct a final choice for the "maximum" coordinate using everything we figured out:
ymaxfinalcoord = [yminusmaxX, yminusmaxY]
Make a final figure of all our work: Mostly taken from: How to fit a 2D ellipse to given points
plt.figure()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('elipse fit to data')
# Plot the original data
plt.scatter(X, Y, linewidths=2, color='black' ,label='Data Points')
# Plot the least squares ellipse
x_coord = np.linspace(xmin-xrng/10,xmax xrng/10,300)
y_coord = np.linspace(ymin-yrng/10,ymax yrng/10,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 x[1] * X_coord * Y_coord x[2] * Y_coord**2 x[3] * X_coord x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2, label='Ellipse Fit')
# Plot the - parabolic solutions of the elipse:
y_parab_soln_plus = yplus( np.linspace( xmin, xmax, 300) )
plt.scatter( np.linspace( xmin, xmax, 300), y_parab_soln_plus, linewidths=2, label='Parab ')
y_parab_soln_minus = yminus( np.linspace( xmin, xmax, 300) )
plt.scatter( np.linspace( xmin, xmax, 300), y_parab_soln_minus, linewidths=2, label='Parab -')
# Plot the maximum value of each parabolic solution:
plt.scatter(yplusmaxX, yplusmaxY, linewidths=2, marker='x', label='yplusmax')
plt.scatter(yminusmaxX, yminusmaxY, linewidths=2, marker='x', label='yminusmax')
plt.scatter(ymaxfinalcoord[0], ymaxfinalcoord[1], s=1360 , marker='x', label='ymaxfinalcoord')
# Plot the naive solution:
plt.scatter(ymaxnaivecoord[0], ymaxnaivecoord[1], linewidths=2, marker='x', label='ymaxnaivecoord')
# Add a legend to the final figure
plt.legend()
#Show all figures:
plt.show()
print('Done')
EDIT: Discussion of the validity of the solution.
This will only work for cases where you have elipses on your image. It should still work if the elipse is cut off (left or right). It should still work if the elipse is rotated. It will work poorly if there are outlier black pixels.
It will not work at all if there is noise in the image, for cases where the elipse is not truely black pixels, and the background is nto truely white. A quick fix to deal with noise would be least squares adding weights. Weights could be used in the least squares fit, to accommodate the darkness of every single pixel above some darkness threshold (could cut most pixels away on greyscale less than some # standard deviations) in the entire image, before performing the fit.
CodePudding user response：
For simplicity I used the matplotlib. However, the code should work with cv2 as well . Maybe these few lines will help you - it was a fun question though :)
import matplotlib.pyplot as plt
import numpy as np
import os
image = plt.imread("lines.png", format="png")
cutoff = .5
rsizex, rsizey = 3,3
color = np.array([1,0,0,])
# find all black pixel by boolean mask
black = ( np.sum( image, axis=2) / 3 ) < cutoff
# x coordinates
# both lines are equivalent, the later is more natural,
# while the first might be easier to understand)
# xcoord = np.tensordot( np.arange( image.shape[0]), np.ones(image.shape[1]), axes= 0)
xcoord = np.reshape( np.arange(image.shape[0]),(-1,1)) * np.ones((1,image.shape[1]) )
# y coordinates
# ycoord = np.tensordot( np.ones( image.shape[0]), np.arange(image.shape[1]), axes= 0)
ycoord = np.reshape( np.arange(image.shape[1]),(1,-1)) * np.ones((image.shape[0],1) )
# coordinates of all black pixel
coord = np.stack( (xcoord, ycoord), axis= 2)[black,:]
# this will find you any pixel with smallest x coordinate
# distance = coord[:,0]
# this loss will find you the pixel with smallest x coordinate and under these
# the pixel which is closest to the center in y.
# Note the second term is always < 1 while coord contains integers
distance = coord[:,0] ** 2 ( ( coord[:,1] - image.shape[1]/2. ) / image.shape[1] ) ** 2
# the resulting coordinate
x,y = coord[ np.argmin(distance), :]
# both blocks are equivalent too, however the second might be easier to understand too,
# while broadcasting is much cleaner !
image[
int( x - rsizex ) : int( x rsizex ),
int( y - rsizey ) : int( y rsizey ),
:,
] = np.reshape( color , (1,1,3) )
# for i, ci in enumerate( [ 1, 0, 0]):
# # paint the block around the found pixel
# image[
# int( x - rsizex ) : int( x rsizex ),
# int( y - rsizey ) : int( y rsizey ),
# i,
# ] = ci
plt.imshow( image )