import cv2 import numpy as np import pandas as pd import matplotlib.pyplot as plt vid = cv2.VideoCapture("GIMBAL.wmv") count = 0 h0 = 29 w0 = 29 h1 = 60 w1 = 60 kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (7, 7)) def I_from_x(x): return (np.identity(2) * (x.T @ x) - (x@x.T)) glangles = [] times = [] frames = [] while True: vid.grab() retval, img = vid.retrieve() if not retval: break # if count < 808: # count += 1 # continue # scale to 01, convert to grayscale image = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) image = image.astype(np.float32) / 255 h, w = (h1, w1) if count >= 800 else (h0, w0) # remove black-hot if count >= 372: image = 1 - image image = (image - np.mean(image)) / np.std(image) center = np.array(image.shape) // 2 image = image[center[0] - h//2:center[0] + h//2, center[1] - w//2:center[1] + w//2] mask = np.zeros(image.shape) mask[image > 0] = 1 mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel) mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, kernel) imageb = image.copy() # standardize again, with respect to the region in the mask only imageb = (imageb - np.mean(imageb)) / np.std(imageb) + 1 imageb[mask == 0] = 0 # center of mass com = np.sum([np.array( idx) * x for idx, x in np.ndenumerate(imageb)], axis=0) / np.sum(imageb) # moment of inertia I = np.sum([I_from_x((np.array(idx) - com).reshape(-1, 1)) * x for idx, x in np.ndenumerate(imageb)], axis=0) / np.sum(imageb) * (h0 * w0) Iw, Iv = np.linalg.eig(I) i1 = np.argmax(Iw) l1 = Iw[i1] v1 = Iv[i1] i2 = np.argmin(Iw) l2 = Iw[i2] v2 = Iv[i2] if (np.cross(v1, v2) < 0): v2 = -v2 glangles.append(np.arctan(-v1[1] / v1[0])) frames.append(count) times.append(count / 29.97) image[mask == 0] = 0 #image = cv2.cvtColor(image, cv2.COLOR_GRAY2RGB) image = img com[0] += center[0] - h//2 com[1] += center[1] - w//2 if count >= 800: scale = 0.6 elif count >= 372: scale = 1.2 else: scale = 1 for v, l, c in zip([v1, v2], [l1, l2], [(0, 255, 0), (255, 0, 0)]): for a in [+scale, -scale]: cv2.line(image, (int(com[1]), int(com[0])), (int( com[1] + a*v[0] * np.sqrt(l) * h / 500), int(com[0] + a*v[1] * np.sqrt(l) * w/500)), c, thickness=1) com = np.round(com).astype(int) for offset in [(0, 0), (-1, 0), (1, 0), (0, -1), (0, 1)]: image[com[0] + offset[0], com[1] + offset[1], :] = (0, 0, 255) # zoom in image = cv2.resize(image, None, fx=2, fy=2, interpolation=cv2.INTER_NEAREST) cv2.imshow("Test", image) cv2.waitKey(1) count += 1 bank = pd.read_csv("bank.csv") bangles = np.array(list(bank["bank"])+[bank["bank"].iloc[-1]],) plt.plot(times, np.rad2deg(glangles)) #plt.plot(times, bangles) plt.show() df = pd.DataFrame.from_dict({ "time": times, "glare_angle": np.rad2deg(glangles), "glare_angle_wrt_horizon": np.rad2deg(glangles) + bangles, "bank_angle": bangles }) # plt.plot(times, np.rad2deg(glangles) + bangles) # plt.show() df["glare_angle_wrt_horizon"] = (np.rad2deg(glangles) + bangles) df["bank_angle"] = bangles df.to_csv("gimbal.csv", index=None)