Полагаю, на вашем месте я хотел бы, чтобы функция приближала края к 90 градусам ... таким образом, чтобы я мог получить период функции sin.Это было бы многошаговым решением.
1) Захватить область, где полосы самые яркие
2) получить значения изображения вдоль линии, которая нормальна к полосам
3) рассчитать кривую наилучшего соответствия.
import cv2
import numpy as np
from matplotlib import pyplot as plt
import math
#function to rotate an image around a point
def rotateImage(image, angle):
image_center = tuple(np.array(image.shape[1::-1]) / 2)
rot_mat = cv2.getRotationMatrix2D(image_center, angle, 1.0)
result = cv2.warpAffine(image, rot_mat, image.shape[1::-1], flags=cv2.INTER_LINEAR)
return result
img = cv2.imread('sin.png', cv2.IMREAD_GRAYSCALE)
imgb = cv2.GaussianBlur(img, (7,7),0)
plt.imshow(imgb)
plt.show()
#grap bright values
v90 = np.percentile(imgb, 90)
#get mask for bright values
msk = imgb >= v90
msk = msk.astype(np.uint8)
plt.imshow(msk)
plt.show()
#get contours
cnts = cv2.findContours(msk, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[-2]
cnt = sorted(cnts, key=cv2.contourArea)
c = np.vstack((cnt[-2], cnt[-1])) #combine 2 largest
((cx, cy), radius) = cv2.minEnclosingCircle(c)
cv2.circle(img, (int(cx), int(cy)), int(radius), 255, 2)
plt.imshow(img)
plt.show()
#circle is just for show...use if you want.
imc = img.copy()
#take top 5 contours(4 would work too ymmv)
angles = []
for i in range(5):
c = cnt[-(i+1)]
ellipse = cv2.fitEllipse(c)
(x, y), (MA, ma), angle = cv2.fitEllipse(c)
cv2.ellipse(imc, ((x, y), (MA, ma), angle), 255, 2)
angles.append(angle)
#average angle of ellipse.
mangle = np.mean(angles)
#goal is create a line normal to the average angle of the ellipses.
#the 0.6 factor here is to just grab the inner region...which will avoid the rapid fall-off of the envelop gaussian-like function.
pt1 = (int(cx + .6*radius*math.cos(math.radians(mangle))), int(cy + .6*radius*math.sin(math.radians(mangle))))
pt2 = (int(cx - .6*radius*math.cos(math.radians(mangle))), int(cy - .6*radius*math.sin(math.radians(mangle))))
#show line
cv2.line(imc, pt1, pt2, 255, 2)
#put fat line on mask...will use this to sample from original image later
im4mask = np.zeros(imc.shape).astype(np.uint8)
cv2.line(im4mask, pt1, pt2, 255, 9)
plt.imshow(imc)
plt.show()
plt.imshow(im4mask)
plt.show()
#now do some rotating(to make the numpy slicing later easier)
imnew = rotateImage(imc, mangle)
plt.imshow(imnew)
plt.show()
im4maskrot = rotateImage(im4mask, mangle)
im4maskrot[im4maskrot > 20] = 255
plt.imshow(im4maskrot)
plt.show()
imgbrot = rotateImage(imgb, mangle)
plt.imshow(imgbrot)
plt.show()
#gather values from original
ys, xs = np.where(im4maskrot == 255)
minx = np.min(xs)
miny = np.min(ys)
maxx = np.max(xs)
maxy = np.max(ys)
print 'x ', minx, maxx
print 'y ', miny, maxy
crop = imgbrot[miny:maxy, minx:maxx]
print crop.shape
plt.imshow(crop)
plt.show()
plt.plot(range(crop.shape[1]), np.mean(crop, axis=0))
#now time to fit a curve.
#first with a gaussian
from scipy import optimize
def test_func(x, a, b, A, mu, sigma):
return A*np.exp(-(x-mu)**2/(2.*sigma**2)) + a * np.sin(b * x)
params, params_covariance = optimize.curve_fit(test_func, np.arange(crop.shape[1]), np.mean(crop, axis=0), p0=[10, 1/15., 60, 2, 150], maxfev=200000000)
print(params)
plt.figure(figsize=(6, 4))
plt.scatter(range(crop.shape[1]), np.mean(crop, axis=0), label='Data')
plt.plot(np.arange(crop.shape[1]), test_func(np.arange(crop.shape[1]), params[0], params[1], params[2], params[3], params[4]), label='Fitted function')
plt.legend(loc='best')
plt.show()
#and without a gaussian...the result is close because of only grabbing a short region.
def test_func(x, a, b, n):
return n + a * np.sin(b * x)
params, params_covariance = optimize.curve_fit(test_func, np.arange(crop.shape[1]), np.mean(crop, axis=0), p0=[10, 1/15., 60], maxfev=200000000)
print(params)
plt.figure(figsize=(6, 4))
plt.scatter(range(crop.shape[1]), np.mean(crop, axis=0), label='Data')
plt.plot(np.arange(crop.shape[1]), test_func(np.arange(crop.shape[1]), params[0], params[1], params[2]), label='Fitted function')
plt.legend(loc='best')
plt.show()
Обратите внимание, что параметр b связан с периодичностью, и оба значения достаточно близки друг к другу (0,0644 и0,0637).Зная это, я бы выбрал более простую подгонку кривой, так как начальные параметры меньше.