Numpy array highest point selection

Question

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.

Lines image

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()

Answer 1

This is a very simple solution if your images are really as you say:

#!/usr/bin/env python3

import cv2
import numpy as np

# Step 1: Load image as greyscale and invert so black is white
im = ~ cv2.imread('5DOwq.png', cv2.IMREAD_GRAYSCALE)

# Step 2: Sum across the rows left to right - black rows will sum to zero...
rowSums = np.sum(im, axis=1)

# Step 3: Find first row in image that doesn't sum to zero, i.e. first row with white pixels
row = np.nonzero(rowSums)[0][0]

# Get middle of all non-zero elements in top row
col = int(np.nonzero(im[row])[0].mean())

print(row,col)

128,283

Graphically, this is what we are doing:

---

Visualisation of the results at each step of the way...

Step 1:

Step 2: rowSums

array([    0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0, 18360, 32640, 41565, 49215, 55845, 43095, 33660,
       29580, 27030, 24225, 21930, 20655, 19380, 17850, 17340, 15810,
       15300, 14790, 13770, 13260, 12495, 12240, 11730, 11220, 10710,
       10200,  9690,  9690,  9180,  8670,  8670,  7905,  7650,  7650,
        7140,  7140,  6885,  6375,  5355,  4845,  4080,  3315,  2550,
        2550,  2040,  1785,  1275,   765,   510,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,   510,   765,  1275,  1530,  2040,  2550,  2550,  3315,
        3825,  4845,  5355,  6375,  6630,  7140,  7140,  7650,  7650,
        8160,  8160,  8670,  9180,  9435, 10200, 10200, 10710, 11220,
       11730, 12240, 12750, 12750, 13770, 14280, 15300, 15810, 16830,
       17850, 18870, 20910, 22440, 23970, 26520, 29580, 33915, 43350,
       56355, 49725, 42075, 32895, 18615,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0,     0,     0,     0,     0,     0,     0,
           0,     0,     0], dtype=uint64)

Step 3: np.nonzero(rowSums)[0]

array([128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140,
       141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153,
       154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166,
       167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 208, 209, 210,
       211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,
       224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236,
       237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249,
       250, 251, 252, 253, 254, 255, 256])

Answer 2

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.

Answer 3

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 =  abs( coord[:,0] ) + abs( coord[:,1] - image.shape[1]/2. ) / image.shape[1] )

# 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 )

However, this only recovers the highest point ( as you stated ) in a naive matter. If you would like to recover the ellipsis as a hole you should take a look at the so called Houghtransform, which is e.g implemented in scikit-learn ( hough_ellipse ).

Fitting the parameters via regession, like in the other answer, is an other option too, but one I would personally not recommend. For example the loss etc can be hard to tune and sometimes erratic for different shapes.

If you need the program working in a more elaborated way and with different shapes you should - of course - use ML, like a CNN. But then we would need a more complete description what you want to archive, where these lines come from and if there is a way to generate data etc.