From fdf95b1734dff70d435fb9e76219ffaf3dc07b85 Mon Sep 17 00:00:00 2001 From: Alex Liberzon Date: Wed, 20 Mar 2024 22:44:49 +0200 Subject: [PATCH 1/2] updated path to tests/test_cavity --- tests/test_pyptv_batch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_pyptv_batch.py b/tests/test_pyptv_batch.py index d78a928..ee258a2 100644 --- a/tests/test_pyptv_batch.py +++ b/tests/test_pyptv_batch.py @@ -3,4 +3,4 @@ def test_pyptv_batch(): # assert cli.cli() == 'CLI template' - pyptv_batch.main('./test_cavity', 10000, 10004) + pyptv_batch.main('./tests/test_cavity', 10000, 10004) From 2c48418f9d50622458a2e82cbff952e824fb78a9 Mon Sep 17 00:00:00 2001 From: Alex Liberzon Date: Wed, 27 Mar 2024 20:51:33 +0200 Subject: [PATCH 2/2] some old graphics file to show differences of rt_is.123456789 --- pyptv/test_calibration.py | 237 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 237 insertions(+) create mode 100755 pyptv/test_calibration.py diff --git a/pyptv/test_calibration.py b/pyptv/test_calibration.py new file mode 100755 index 0000000..312f6f5 --- /dev/null +++ b/pyptv/test_calibration.py @@ -0,0 +1,237 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +""" +Created on Fri Nov 17 09:57:49 2017 + +@author: ron + + +Here is a script that is meant to help in the task of evaluating the quality of +PyPTV calibration. + +Once a PyPTV experiment folder is ready and a calibration is established, the +evaluation here is made by comparing known points of calibration (i.e. calblock) +points, with points that were determined using images of the calibration target +(i.e. dt_lsq points). To generate the dt_lsq points load the calibration images +as the ones to analyze first. Then process the images with: + + image coords -> corespondeces -> 3D Positions + +The script here is used by loading the point files with the functions: +read_dt_lsq(), and read_calblock(). After that use the function pair_cal_points() +to match points from both sets. + +The evaluation itself is made by first plotting the points in 3D with +plot_cal_points(). The distribution of errors can the be examined with the +plot_cal_err_histogram() function +""" + +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + + + + +def read_dt_lsq(file_path): + """ + will read a PyPTV dt_lsq file and return the points as a list of + numpy arrays + + inputs + ====== + file_path (string) - absolute path to dt_lsq file + + output + ====== + points (list) - list of numpy (3,1) arrays with (x,y,z) coordinates + """ + f = open(file_path,'r') + N_particles = int(f.readline().strip()) + points = [] + + for i in range(N_particles): + l = f.readline().strip().split() + point = np.array([l[1], l[2], l[3]], dtype=float) + points.append(point) + + f.close() + + return points + + + + +def read_calblock(file_path): + """ + will read a PyPTV calbloack file and return the points as a list of + numpy arrays + + inputs + ====== + file_path (string) - absolute path to dt_lsq file + + output + ====== + points (list) - list of numpy (3,1) arrays with (x,y,z) coordinates + """ + f = open(file_path,'r') + a = f.readlines() + f.close() + points = [] + + for i in range(len(a)): + l = a[i].strip().split() + try: + point = np.array([l[1], l[2], l[3]], dtype=float) + except: + print('last data', l) + raise ValueError('bad line in calblock file') + points.append(point) + + return points + + + + +def pair_cal_points(calblock_pnts, dt_lsq_pnts, max_dist = 3.0): + ''' + will determine pairs of points from the dt_lsq file and the known calblock + file. for each point in the dt_lsq file, will find the closest point to it + from the calblock points. + + inputs + ====== + calblock_pnts (list) - a list of array(3,1) points from a calblock file + dt_lsq_pnts (list) - a list of array(3,1) points for a dt_lsq file + max_dist (float) - the maximum distance that can be regarded a pair + + output + ====== + pairs_list (list) - a list of pairs of points. the first is a calbclock + point and the second a dt_lsq point + + ''' + N_cb = len(calblock_pnts) + N_dt = len(dt_lsq_pnts) + N_pairs = min(N_cb, N_dt) + + dist_mat = np.zeros( (N_cb, N_dt) ) + index_mat = np.zeros( (N_cb, N_dt), dtype=[ ('i', 'i4'),('j','i4' )]) + for i in range(dist_mat.shape[0]): + for j in range(dist_mat.shape[1]): + dist_mat[i,j] = np.linalg.norm(calblock_pnts[i] - dt_lsq_pnts[j]) + index_mat[i,j] = (i,j) + + pairs_list = [] + for i in range(N_pairs): + d = np.amin(dist_mat) + if d < max_dist: + w = np.where(dist_mat == np.amin(dist_mat)) + i_ = index_mat['i'][w[0][0], w[1][0]] + j_ = index_mat['j'][w[0][0], w[1][0]] + pairs_list.append( (calblock_pnts[i_], + dt_lsq_pnts[j_]) ) + + dist_mat = np.delete(dist_mat, w[0][0], axis=0) + dist_mat = np.delete(dist_mat, w[1][0], axis=1) + index_mat = np.delete(index_mat, w[0][0], axis=0) + index_mat = np.delete(index_mat, w[1][0], axis=1) + else: break + return pairs_list + + + + + +def plot_cal_points(pairs_list): + ''' + plot a 3D scatter plot of calblock points (red) and dt_lsq points (blue). + + input + ===== + pairs_list (list) - output from pair_cal_points() + + output + ====== + fig, ax - matplotlib figure and axis objets + ''' + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + for p in pairs_list: + ax.plot([p[0][0]], [p[0][2]], [p[0][1]], 'xr') + ax.plot([p[1][0]], [p[1][2]], [p[1][1]], 'xb') + + ax.set_xlabel('X') + ax.set_ylabel('Z') + ax.set_zlabel('Y') + + return fig, ax + + + + +def plot_cal_err_histogram(pairs_list): + ''' + plot a 3D scatter plot of calblock points (red) and dt_lsq points (blue). + + input + ===== + pairs_list (list) - output from pair_cal_points() + + output + ====== + fig, ax - matplotlib figure and axis objets + ''' + + dx,dy,dz = [],[],[] + + for p in pairs_list: + dx.append(p[0][0] - p[1][0]) + dy.append(p[0][1] - p[1][1]) + dz.append(p[0][2] - p[1][2]) + + fig,ax = plt.subplots() + + lbls = [r'x',r'y',r'z'] + for e,lst in enumerate([dx,dy,dz]): + m,s = np.mean(lst), np.std(lst) + h=ax.hist(lst,bins=8,histtype= 'step', lw=3, + label=r'$\langle %s \rangle=%0.3f, \sigma_{%s}=%0.3f$'%(lbls[e],m,lbls[e],s)) + #h = np.histogram(lst,bins=10) + #x,y = (h[1][:-1] + h[1][1:])*0.5 , h[0] + #ax.plot(x,y, '-o', lw=2, + #label=r'$%s:$ $\mu=%0.0e, \sigma=%0.1e$'%(lbls[e],m,s)) + ax.legend(loc='best') + return fig, ax + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file