From 830ed04ff886c083afc42a9b19a06e988fd5d7bd Mon Sep 17 00:00:00 2001 From: Ron Shnapp Date: Sun, 14 Jan 2024 00:38:35 +0200 Subject: [PATCH] extended zolof model --- myptv/TsaiModel/camera.py | 33 +--- myptv/extendedZolof/__init__.py | 12 ++ myptv/extendedZolof/calibrate.py | 103 ++++++++++++ myptv/extendedZolof/camera.py | 280 +++++++++++++++++++++++++++++++ myptv/utils.py | 137 +++++++++------ setup.py | 2 +- 6 files changed, 482 insertions(+), 85 deletions(-) create mode 100644 myptv/extendedZolof/__init__.py create mode 100644 myptv/extendedZolof/calibrate.py create mode 100755 myptv/extendedZolof/camera.py diff --git a/myptv/TsaiModel/camera.py b/myptv/TsaiModel/camera.py index c4f5e42..556c44d 100755 --- a/myptv/TsaiModel/camera.py +++ b/myptv/TsaiModel/camera.py @@ -43,6 +43,7 @@ """ +from myptv.utils import Cal_image_coord import os from math import sin, cos @@ -362,38 +363,6 @@ def plot_3D_epipolar_line(self, eta, zeta, zlims, ax=None, color=None): -class Cal_image_coord(object): - ''' - A class used for reading the calibration image files. This is called - by the camera class if given a filename with calibration points. - ''' - - def __init__(self, fname): - ''' - input - - fname - String, the path to your calibration point file. The file is - holds tab separated values with the meaning of: - [x_image, y_image, x_lab, y_lab, z_lab] - ''' - self.image_coords = [] - self.lab_coords = [] - self.fname = fname - self.read_file() - - - def read_file(self): - - with open(self.fname) as f: - lines = f.readlines() - self.N_points = len(lines) - - for ln in lines: - ln_ = ln.split() - self.image_coords.append([float(ln_[0]), float(ln_[1])]) - self.lab_coords.append( [float(ln_[2]), - float(ln_[3]), - float(ln_[4])] ) - f.close() diff --git a/myptv/extendedZolof/__init__.py b/myptv/extendedZolof/__init__.py new file mode 100644 index 0000000..4c0f6d2 --- /dev/null +++ b/myptv/extendedZolof/__init__.py @@ -0,0 +1,12 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun March 20 2020 + +@author: ron + +MyPTV init file +""" + + + + \ No newline at end of file diff --git a/myptv/extendedZolof/calibrate.py b/myptv/extendedZolof/calibrate.py new file mode 100644 index 0000000..e0a383a --- /dev/null +++ b/myptv/extendedZolof/calibrate.py @@ -0,0 +1,103 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Nov 1 12:15:21 2019 + +@author: ron + + +calibration module for Extended Zolof camera instances. This is used to +obtain the [A], [B], and O Extended Zolof model parameters. + +""" + +from numpy import array +from numpy import sum as npsum +from numpy.linalg import lstsq, norm +from scipy.optimize import minimize + +from myptv.extendedZolof.camera import camera_extendedZolof +from myptv.utils import line, get_nearest_line_crossing + + + + + +class calibrate_extendedZolof(camera_extendedZolof): + ''' + This object is used to calibrate cameras against a given list + of lab and camera point coordinates. + ''' + + def __init__(self, camera, x_list, X_list): + ''' + Given a list of 2D points, x=(x,y), and a list of 3D point X=(X,Y,Z), + we assume that given a point X, we can compute x by a polynomial + of degree 3, as - + + x = A0 + A1*X + A2*Y + A3*Z + + A4*X^2 + A5*Y^2 + A6*Z^2 + A7*XY + A8*YZ + A9*ZX + A108XYZ + A11*XY^2 + A12*XZ^2 + A13*YX^2 + A14*YZ^2 + A15*ZX^2 + A16*ZY^2 + ''' + self.cam = camera + self.A = [[0.0 for i in range(17)] for j in [0,1]] + self.B = [[0.0 for i in range(10)] for j in [0, 1, 2]] + self.x_list = x_list + self.X_list = X_list + + + + + def calibrate(self): + ''' + Given a list of points, x and X, this function attempts to determine + the A coefficients. + ''' + # 1) finding the A coefficients - + XColumns = [self.cam.get_XCol(Xi) for Xi in self.X_list] + res = lstsq(XColumns, self.x_list, rcond=None) + self.A = res[0] + + # 2) finding the best camera center - + line_list = [] + for i in range(0, len(self.X_list)): + O, e = self.get_ray_from_x(self.x_list[i], X0=self.X_list[i]) + line_list.append(line(O, e)) + self.O = get_nearest_line_crossing(line_list) + + # 3) finding the unit vector for each X - + r_list = [] + for Xi in self.X_list: + r = (Xi - self.O)/norm(Xi - self.O) + r_list.append(r) + + # 4) finding the B coefficients - + xColumns = [self.cam.get_xCol(xi) for xi in self.x_list] + res = lstsq(xColumns, r_list, rcond=None) + self.B = res[0] + + self.cam.O = self.O + self.cam.A = self.A + self.cam.B = self.B + + + + def get_ray_from_x(self, x, X0=None): + ''' + Given a point in 2D image space, this function returns a line in 3D + that passes through this point. The line is represented with six + parameters: one point in 3D, O, and one unit vector in 3D, e. + ''' + + func = lambda X: sum((array(self.projection(X)) - array(x))**2) + + if X0 is None: + X0 = array([0,0,0]) + + X02 = array(X0) + array([1,1,1]) + + O = minimize(func, X0).x + dX = minimize(func, X02).x + e = (O-dX)/sum((O-dX)**2)**0.5 + + return O, e + \ No newline at end of file diff --git a/myptv/extendedZolof/camera.py b/myptv/extendedZolof/camera.py new file mode 100755 index 0000000..00d91e2 --- /dev/null +++ b/myptv/extendedZolof/camera.py @@ -0,0 +1,280 @@ +# -*- coding: utf-8 -*- +#!/usr/bin/env python3 + +""" +Created on Jan 13 2024 + +@author: ron + + +containts the camera class that represents the Extended Zolof 3D camera model. +The class handles the transformation from camera space coordinates to +lab coordinates using the extended Zolof model, with a third-order back-and- +forth polynomial transformation. + + + +The original forward transformation model: + +The Extended Zolof 3D model is an extension of the existing Zolof polynomial +method of 3D transformation. In the original Zolof model, each 3D point in +lab-space, X, is related to a 2D camera-space point, x, through a third degree +poynomial. Specifically, let X = (X1, X2, X3) and x = (eta, zeta). The +transformation from X to x is: + + x = [A] * P(X); + +In MyPTV, P(X) is the following third order polynomial with 17 terms: + + P(X) = 1 + X1 + X2 + X3 + X1^2 + X2^2 + X3^2 + X1*X2 + X2*X3 + X3*X1 + + X1*X2^2 + X1*X3^2 + X2*X1^2 + X2*X3^2 + X3*X1^2 + X3*X2^2 + X1*X2*X3 + +and [A] is a vector of 17 coefficients: + + [A] = A1, A2, ..., A17 . + +Note that unlike other models, the Zolof model does not assume anything about +the physics of the camera, and simply applied a least-squares fitting to find +the coefficients [A] that reduce the transformaiton error. The polynomial used +in MyPTV was tested to give good results for standard PTV experiments. + + + +The extension to backward transformation: + +The original Zolof model can only perform a lab-to-camera transformation. +However, to complete a PTV analysis, we need to be able to transform camera +space coordinates back into lab space coordinates, and this does not exist in +the original Zolof model. Therefore, in MyPTV we extended this model to +include also the inverse, camera-to-lab transformation. + +To do so we use an approach that combines the pin-hole model with the Zolof +polynomial. Specifically, we assume that each point in camera-space (e.g. +each camera pixel) corresponds to a stright physical line in lab-space; these +lines are called epipolar lines. Furthermore, we assume that all the epipolar +lines for a certain camera cross each other in a single point. These two +assumptions were tested with two separate experiments, which gave good results. +With the two assumptions given above, it is possible to define the backwards, +namely, camera-to-lab transformation as a function that given coordinates +in camera space, returns an epipolar line. + +With the second assumption above, epipolar lines belonging to a camera can +be defined using a direction vector. Thus, the backwards, camera-to-lab +transformation is made as follows. Let O = (O1, O2, O3) denote the point in +which all epipolar lines cross, and let r = (r1, r2, r3) denote a direction +vector in lab-space. Points laying along an epipolar line can thus be written +as: + + l(x) = O + r(x) * a + +where a is a free parameter running from -infinity to +infinity. The +camera-to-lab transformation is achieved by: + + r(x) = [B] G(x) + +where G(x) is a third degree polynomial with the following 10 terms + + G(x) = 1 + x1 + x2 + x1^2 + x2^2 + x1*x2 + x1^2*x2 + x2^2*x1 + x1^3 + x2^3 + +and where [B] is a vector of 10 polynomial coeficients: + + [B] = B1, B2, ... , B10 . + + + + +The advantage of the Extended Zolof model: + +The advantage of the Extended Zolof model is that 1) it is fully defined and +and perfectly fits to be used in PTV experiments because it has the forward +and backwards transformations; 2) it is straightforward to determine the +model parameters in the calibration process. Specifically, the problem can be +formulated as a normal least-squares minimization, which gives out fast and +unique results without resorting to numerial minimization. This is true for +all the model parameters, [A], [B], and O. + +""" + +from myptv.utils import Cal_image_coord + +from numpy import dot, array, zeros +import os + + + + + + +class camera_extendedZolof(object): + ''' + an object that holds the calibration information for each camera using the + Extended Zolof model. + It can be used to: + 1) obtain image coordinates from given lab coordinates. + 2) vice versa if used together with other cameras at + other locations (img_system.stereo_match). + + input: + name - string name for the camera + resolution - tuple (2) two integers for the camera pixels + cal_points_fname - path to a file with calibration coordinates for the cam + ''' + + def __init__(self, name, cal_points_fname = None): + ''' + inputs: + name - (string) name of the camera and name of the file in which + it is stored on the disk. + + cal_points_fname (string) - path of a file that stores the calibration + points data. Optional. + ''' + + self.name = name + self.O = zeros(3) + 1. # camera location + self.A = array([[0.0 for i in range(17)] for j in [0,1]]).T + self.B = array([[0.0 for i in range(10)] for j in [0, 1, 2]]).T + + if cal_points_fname is not None: + cic = Cal_image_coord(cal_points_fname) + self.image_points = cic.image_coords + self.lab_points = cic.lab_coords + + + + def __repr__(self): + + ret = ('Extended Zolof model camera instace. ' + + self.name + '\n O: ' + str(self.O)) + return ret + + + + def get_XCol(self, X): + ''' + Given a point in 3D lab-space, this method returns its 17 P(X) + polynomial terms. + ''' + X1,X2,X3 = X[0],X[1],X[2] + XColumn = [1.0, X1, X2, X3, + X1**2, X2**2, X3**2, X1*X2, X2*X3, X3*X1, X1*X2*X3, + X1*X2**2, X1*X3**2, X2*X1**2, X2*X3**2, X3*X1**2, X3*X2**2] + return XColumn + + + + def get_xCol(self, x): + ''' + Given a point in 2D camera-space, this method returns its 10 G(x) + polynomial terms + ''' + x1, x2 = x[0], x[1] + xColumn = [1.0, x1, x2, x1**2, x2**2, x1*x2, + x1**2*x2, x2**2*x1, x1**3, x2**3] + return xColumn + + + + def projection(self, X): + ''' + Given a point in 3D, X, this method returns its 2D camera-space + projection. + ''' + XColumn = self.get_XCol(X) + res = dot(XColumn, self.A) + return [res[0], res[1]] + + + + def get_r(self, eta, zeta): + ''' + Given a point in 2D, x, this method returns its 3D direction vector + ''' + x = [eta, zeta] + xColumn = self.get_xCol(x) + res = dot(xColumn, self.B) + return [res[0], res[1], res[2]] + + + + def get_r_ori(): + msg1 = 'The extended Zolof model is not yet capable of returning' + msg2 = 'particle orientation results. Use the Tsai model instead.' + raise TypeError(msg1 + msg2) + + + + def save(self, dir_path = ''): + full_path = os.path.join(dir_path, self.name) + + f = open(full_path, 'w') + f.write('extendedZolof model camera\n') + f.write(self.name+'\n') + + S = '' + for i in range(len(self.O)): + S += 'O %s \t %s \n'%(i, self.O[i]) + + for i in range(len(self.A)): + for j in range(len(self.A[i])): + S += 'A %s %s \t %s \n'%(i, j, self.A[i][j]) + + for i in range(len(self.B)): + for j in range(len(self.B[i])): + S += 'B %s %s \t %s \n'%(i, j, self.B[i][j]) + + f.write(S) + f.close() + + + + def load(self, dir_path): + ''' + will load camera data from the hard disk + ''' + full_path = os.path.join(dir_path, self.name) + + f = open(full_path, 'r') + model = f.readline() + modelName = model.split()[0] + if modelName != 'extendedZolof': + msg = 'Camera file is not an extendedZolof camera (%s)'%modelName + raise ValueError(msg) + name = f.readline() + + lines = f.readlines() + f.close() + + for i in range(len(lines)): + lines[i] = lines[i].strip().split() + + O_vals = list(filter(lambda l: l[0]=='O', lines)) + for v in O_vals: + ind = int(v[1]) + self.O[ind] = float(v[-1]) + + A_vals = list(filter(lambda l: l[0]=='A', lines)) + Ashape = (int(max(A_vals, key=lambda x: float(x[1]))[1])+1, + int(max(A_vals, key=lambda x: float(x[2]))[2])+1) + self.A = array([[0.0 for i in range(Ashape[0])] for j in range(Ashape[1])]).T + for v in A_vals: + i, j = int(v[1]), int(v[2]) + self.A[i][j] = float(v[-1]) + + B_vals = list(filter(lambda l: l[0]=='B', lines)) + Bshape = (int(max(B_vals, key=lambda x: float(x[1]))[1])+1, + int(max(B_vals, key=lambda x: float(x[2]))[2])+1) + self.B = array([[0.0 for i in range(Bshape[0])] for j in range(Bshape[1])]).T + for v in B_vals: + i, j = int(v[1]), int(v[2]) + self.B[i][j] = float(v[-1]) + + + + + + + + + + diff --git a/myptv/utils.py b/myptv/utils.py index ab55810..8f28fc4 100644 --- a/myptv/utils.py +++ b/myptv/utils.py @@ -83,6 +83,50 @@ def point_line_dist(O,r,P): +class line(object): + + def __init__(self, O, r): + ''' + A helper class for calculating line and point distances. + O - a point through which the line passes + r - a vector of the line direction + ''' + self.O = O + self.r = r + + def distance_to_point(self, P): + ''' + Returns the least distance between the line and a given point P. + ''' + s1 = sum([self.r[i]*(P[i]-self.O[i]) for i in [0,1,2]]) + s2 = sum([ri**2 for ri in self.r]) + amin = s1 / s2 + + dmin = sum([(self.O[i]-P[i]+amin*self.r[i])**2 for i in [0,1,2]])**0.5 + + return dmin, amin, self.O + amin*self.r + + + + +def get_nearest_line_crossing(line_list): + ''' + Given a list of line objects, this function find a point that minimizes + the sum of the distances to it. + ''' + from scipy.optimize import minimize + func = lambda P: sum([l.distance_to_point(P)[0] for l in line_list]) + P = minimize(func, array([0,0,0])).x + return P + + + + + + + + + def fit_polynomial(x, y, n): ''' A polynomial of degree n is written as: @@ -119,6 +163,47 @@ def fit_polynomial(x, y, n): + + + +class Cal_image_coord(object): + ''' + A class used for reading the calibration image files. This is called + by the camera class if given a filename with calibration points. + ''' + + def __init__(self, fname): + ''' + input - + fname - String, the path to your calibration point file. The file is + holds tab separated values with the meaning of: + [x_image, y_image, x_lab, y_lab, z_lab] + ''' + self.image_coords = [] + self.lab_coords = [] + self.fname = fname + self.read_file() + + + def read_file(self): + + with open(self.fname) as f: + lines = f.readlines() + self.N_points = len(lines) + + for ln in lines: + ln_ = ln.split() + self.image_coords.append([float(ln_[0]), float(ln_[1])]) + self.lab_coords.append( [float(ln_[2]), + float(ln_[3]), + float(ln_[4])] ) + f.close() + + + + + + class match_calibration_blobs_and_points(object): ''' This is a tool used to matxh between a calibration "target file" and @@ -209,58 +294,6 @@ def plot_projections(self): -#============================================================================= -#Line intersects. Not in use at the moment - only, testing. - - -def find_point_nearest_to_lines(line_list): - ''' - Will use least squares in order to solve the problem of finding the - point that is the closest to several lines. (Note that for two lines there - is an analytical solution that is used in the function line_dist()). - - This function uses scipy.optimize.least_squares() ! - - input - - line_list - a list of tuples (O, r) where O and r are the origin and - direction vectors (numpy arrays, 3) that define a line in - 3d space. - - returns - - P - the point in 3d space that minimizes the distance to the lines in the - list. - ''' - from scipy.optimize import least_squares - - def dist_sum(P): - dist = 0 - for O,r in line_list: - dist += point_line_dist(O, r, P) - return dist - - x0 = line_dist(line_list[0][0], line_list[0][1], - line_list[1][0], line_list[1][1])[1] - P = least_squares(dist_sum, x0, method='dogbox', ftol=1e-8) - return P - - - - - -def nearest_intersect(lines): - import numpy as np - I = np.eye(3) - - S1 = 0 - S2 = 0 - for i in range(len(lines)): - Oi = np.array([lines[i][0]]) - ri = lines[i][1] - S1 += I - np.dot(Oi.T, Oi) - S2 += np.dot(I - np.dot(Oi.T, Oi), ri) - return np.linalg.lstsq(S1, S2, rcond=None) - -#============================================================================= diff --git a/setup.py b/setup.py index 12b350e..97ae969 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ setup( name='myptv', - packages=find_packages(include=['myptv', 'myptv.fibers', 'myptv.TsaiModel']), + packages=find_packages(include=['myptv', 'myptv.fibers', 'myptv.TsaiModel', 'myptv.extendedZolof']), version='1.0.0', description='A 3D Particle Tracking Velocimetry library', install_requires=['numpy', 'scipy', 'scikit-image','pandas','matplotlib','pyyaml', 'tk', 'Pillow'],