one cam combine

config/cfg_3freq_wafer_1.json

@@ -0,0 +1,16 @@
+{
+    "cam_num":1,
+    "screen_square_size": 29,
+    "screen_chessboard_size": [10, 7],
+    "world_chessboard_size": [8, 8],
+    "d": -20.7,
+    "num_freq": 3,
+    "screenPixelSize": [1920, 1080],
+    "screenorigin_point": [1320, 300],
+    "Axis_conversion": [[-1, 1]],
+    "k": 32.0,
+    "pixelSize_mm": 0.363745,
+    "grid_spacing": 1,
+    "cam_params": "cam_params_1.pkl",
+    "screen_params":"screen_params_1.pkl"
+}

config/cfg_3freq_wafer_4.json

@@ -0,0 +1,17 @@
+{
+    "cam_num":4,
+    "screen_square_size": 25.31,
+    "screen_chessboard_size": [10, 6],
+    "world_chessboard_size": [6, 6],
+    "world_square_size": 35.2,
+    "d": -2.5,
+    "num_freq": 3,
+    "screenPixelSize": [1600, 1200],
+    "screenorigin_point": [649, 550],
+    "Axis_conversion": [[-1, -1]],
+    "k": 32.0,
+    "pixelSize_mm": 0.253125,
+    "grid_spacing":1,
+    "cam_params": "cam_params_4.pkl",
+    "screen_params":"screen_params_4.pkl"
+}

config/cfg_3freq_wafer_matlab.json

@@ -0,0 +1,16 @@
+{
+    "cam_num":1,
+    "screen_square_size": 9,
+    "screen_chessboard_size": [10, 7],
+    "world_chessboard_size": [8, 8],
+    "d": 10.04,  
+    "num_freq": 3,
+    "screenPixelSize": [2160, 1215],
+    "screenorigin_point": [1512.5,121.5],
+    "Axis_conversion": [[-1, 1]],
+    "k": 32.0,
+    "pixelSize_mm": 0.095955,
+    "grid_spacing": 1,
+    "cam_params": "cam_params_1.pkl",
+    "screen_params":"screen_params_1.pkl"
+}

pmd.py

@@ -24,6 +24,9 @@ import cv2
 from src.eval import get_eval_result
 import pickle
 from collections import defaultdict
+from scipy.io import loadmat, savemat
+
+
 
 def pmdstart(config_path, img_folder):
     start_time = time.time()
@@ -42,15 +45,15 @@ def main(config_path, img_folder):
     os.chdir(current_dir)
 
     cfg = json.load(open(config_path, 'r'))
-    n_cam = 4
+    n_cam = cfg['cam_num']
     num_freq = cfg['num_freq']
     save_path = 'debug'
     debug = False
     grid_spacing = cfg['grid_spacing']
     num_freq = cfg['num_freq']
-    smooth = True
-    align = True
-    denoise = True
+    smooth = False
+    align = False
+    denoise = False
     #cammera_img_path = 'D:\\data\\four_cam\\calibrate\\calibrate-1008'
     screen_img_path = 'D:\\data\\four_cam\\calibrate\\cam3-screen-1008'
     cammera_img_path = 'D:\\data\\four_cam\\calibrate\\calibrate-1016'
@@ -60,8 +63,10 @@ def main(config_path, img_folder):
     print("\n1. 相机标定")
     preprocess_start = time.time()
 
-    #cam_para_path = os.path.join('config', cfg['cam_params'])
-    cam_para_path = 'D:\\code\\pmd-python\\config\\cam_params.pkl'
+    cam_para_path = os.path.join(current_dir, 'config', cfg['cam_params'])
+    print('cam_para_path = ', cam_para_path)
+    print('current_dir = ', current_dir)
+    #cam_para_path = 'D:\\code\\pmd-python\\config\\cam_params.pkl'
     if os.path.exists(cam_para_path):
     #if False:
         with open(cam_para_path, 'rb') as pkl_file:
@@ -77,8 +82,8 @@ def main(config_path, img_folder):
             #print('cam_img_path = ', cam_img_path)
             cam_param_raw = calibrate_world(cam_img_path, i, cfg['world_chessboard_size'], cfg['world_square_size'], debug=0)
             cam_params.append(cam_param_raw)
-        with open(cam_para_path, 'wb') as pkl_file:
-           pickle.dump(cam_params, pkl_file) 
+        # with open(cam_para_path, 'wb') as pkl_file:
+        #    pickle.dump(cam_params, pkl_file) 
 
     print("\n2. 屏幕标定")
     screen_cal_start = time.time()
@@ -91,10 +96,10 @@ def main(config_path, img_folder):
             screen_params = pickle.load(pkl_file)[0]
     else:
         screen_params = calibrate_screen(screen_img_path, cam_params[3]['camera_matrix'], cam_params[3]['distortion_coefficients'], cfg['screen_chessboard_size'], cfg['screen_square_size'], debug=0)
-        with open(screen_para_path, 'wb') as pkl_file:
-            pickle.dump([screen_params], pkl_file)
+        # with open(screen_para_path, 'wb') as pkl_file:
+        #     pickle.dump([screen_params], pkl_file)
+    
     
-    screen_to_world = map_screen_to_world(screen_params, cam_params[3])  
     screen_cal_end = time.time()
     print(f"   完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
     print(f"   耗时: {screen_cal_end - screen_cal_start:.2f} 秒")
@@ -105,18 +110,109 @@ def main(config_path, img_folder):
     binary_masks = []
     for cam_id in range(n_cam):
         print('cam_id = ', cam_id)
-        white_path = os.path.join(img_folder, f'{cam_id}_frame_24.bmp') 
-        binary = get_white_mask(white_path, bin_thresh=12, debug=0)
+        #white_path = os.path.join(img_folder, f'{cam_id}_frame_24.bmp') 
+        white_path = 'D:\\code\\code of PMD\\code of PMD\\picture\\white.bmp'
+        if n_cam == 3:
+            binary = get_white_mask(white_path, bin_thresh=12, debug=0)
+            
+        elif n_cam == 1:
+            #angle_rad, binary = find_notch(white_path, n_cam, debug=0)
+            binary = get_white_mask(white_path, bin_thresh=12, debug=0)
+
         binary_masks.append(binary)
-        phases = extract_phase(img_folder, cam_id, binary, cam_params[cam_id]['camera_matrix'], cam_params[cam_id]['distortion_coefficients'], num_freq=num_freq)
-        x_un, y_un = unwrap_phase(phases, save_path, num_freq, debug=0)
-        x_uns.append(x_un)
-        y_uns.append(y_un)
+
+        #phases = extract_phase(img_folder, cam_id, binary, cam_params[cam_id]['camera_matrix'], cam_params[cam_id]['distortion_coefficients'], num_freq=num_freq)
+        #x_un, y_un = unwrap_phase(phases, save_path, num_freq, debug=0)
+        
+        #x_uns.append(x_un)
+        #y_uns.append(y_un)
 
     phase_end = time.time()
     print(f"   完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
     print(f"   耗时: {phase_end - phase_start:.2f} 秒")
- 
+    
+    # matlab align
+    #cam_params = loadmat('D:\\code\\code of PMD\\code of PMD\\calibration\\calibrationSessionworld.mat')
+    #screen_params = loadmat('D:\\code\\code of PMD\\code of PMD\\calibration\\calibrationSessionscreen.mat')
+   
+    #print('cam_params = ', cam_params)
+    #print('screen_params = ', screen_params)
+
+    cam_params = []
+    screen_params = []
+    mtx = np.array([[6944.89018564351, 0, 2709.44699784446], [0, 6942.91962497637, 1882.05677580185], [0,0,1]])
+    R = cv2.Rodrigues(np.array([0.0732148059333282,-0.265812028310130,-0.0532640604086260]).reshape(3,1))[0]
+    T = np.array([-15.179977,-42.247126,246.92182]).reshape(3,1)
+    
+    cam_calibration_data = {
+        'camera_matrix': mtx,
+        'distortion_coefficients': 0,
+        'rotation_matrix': R,
+        'translation_vector': T,
+        'error': 0
+    }
+
+    mtx = np.array([[6944.89018564351, 0, 2709.44699784446], [0, 6942.91962497637, 1882.05677580185], [0,0,1]])
+    R = np.array([[0.96514996,-0.042578806,0.25821037],[0.029285983,0.99805061,0.055111798],[-0.26005361,-0.045629206,0.96451547]])
+    T = np.array([30.1970,-49.3108,507.4424]).reshape(3,1)
+    
+    screen_calibration_data = {
+        'screen_matrix': mtx,
+        'screen_distortion': 0,
+        'screen_rotation_matrix': R,
+        'screen_translation_vector': T,
+        'screen_error': 0
+    }
+
+    
+    cam_params.append(cam_calibration_data)
+    #screen_params.append(screen_calibration_data)
+    screen_params = screen_calibration_data
+
+    x_un = loadmat('D:\\code\\code of PMD\\code of PMD\\phase_information\\x_phase_unwrapped.mat')
+    y_un = loadmat('D:\\code\\code of PMD\\code of PMD\\phase_information\\y_phase_unwrapped.mat')
+
+    x_un = x_un['x_phase_unwrapped']
+    y_un = y_un['y_phase_unwrapped']
+
+    x_uns, y_uns = [], []
+    x_uns.append(x_un)
+    y_uns.append(y_un)
+   
+    #fig, axes = plt.subplots(1, 2, figsize=(12, 6))
+
+
+        # 第一个子图
+    # cax0 = axes[0].imshow(x_un)
+    # axes[0].set_title('x_phase_unwrapped')
+    # axes[0].set_xlabel('X Axis')
+    # axes[0].set_ylabel('Y Axis')
+    # fig.colorbar(cax0, ax=axes[0])
+
+    #     # 第二个子图
+    # cax1 = axes[1].imshow(y_un)
+    # axes[1].set_title('y_phase_unwrapped')
+    # axes[1].set_xlabel('X Axis')
+    # axes[1].set_ylabel('Y Axis')
+    # fig.colorbar(cax0, ax=axes[1])
+
+    #     # 调整子图之间的间距
+    # plt.tight_layout()
+        #plt.savefig(os.path.join(save_path, "phase_unwrapped.png"))
+    #plt.show()
+    
+    #return 0
+    print('screen_params = ', screen_params)
+
+    if n_cam == 1:
+        screen_to_world = map_screen_to_world(screen_params, cam_params[0]) 
+    elif n_cam == 4:
+        screen_to_world = map_screen_to_world(screen_params, cam_params[3])  
+    else:
+        print('camera number should be 1 or 4')
+        return 0
+    
+
     print("\n4. 获得不同坐标系下点的位置")
     get_point_start = time.time()
     total_cloud_point = np.empty((0, 3))
@@ -132,7 +228,7 @@ def main(config_path, img_folder):
                     'u_p': u_p, 'v_p': v_p}
         screen_points = get_screen_points(point_data, x_uns[i], y_uns[i], screen_params, screen_to_world, cfg, save_path, i, debug=debug)
         #plot_coords(world_points, camera_points, screen_points)
-        z_raw, aligned, smoothed, denoised, gradient_xy = reconstruction_cumsum(world_points, camera_points, screen_points, save_path, i, debug=0, smooth=smooth, align=align, denoise=denoise)
+        z_raw,  gradient_xy = reconstruction_cumsum(world_points, camera_points, screen_points, debug=0, smooth=smooth, align=align, denoise=denoise)
         
         z_raw_xy = np.round(z_raw[:, :2]).astype(int)
         
@@ -161,12 +257,13 @@ def main(config_path, img_folder):
         # z_raw_aligned = non_boundary_points @ rotation_matrix.T
         # z_raw_aligned[:,2] = z_raw_aligned[:,2] - np.mean(z_raw_aligned[:, 2])
 
-        z_raw_aligned = non_boundary_points
+        #z_raw_aligned = non_boundary_points
+        z_raw_aligned, _  = align2ref(non_boundary_points)
 
         #non_boundary_points = smoothed
-        write_point_cloud(os.path.join(img_folder, str(i) + '_cloudpoint.txt'), np.round(z_raw_aligned[:, 0]), np.round(z_raw_aligned[:, 1]),  1000*z_raw_aligned[:, 2])
+        write_point_cloud(os.path.join(img_folder, str(i) + '_cloudpoint.txt'), np.round(z_raw_aligned[:, 0]), np.round(z_raw_aligned[:, 1]),  z_raw_aligned[:, 2])
         np.savetxt(os.path.join(img_folder, str(i) + '_gradient.txt'), gradient_xy, fmt='%.10f', delimiter=',')
-        total_cloud_point = np.vstack([total_cloud_point, np.column_stack((z_raw_aligned[:, 0], z_raw_aligned[:, 1], 1000*z_raw_aligned[:, 2]))])
+        total_cloud_point = np.vstack([total_cloud_point, np.column_stack((z_raw_aligned[:, 0], z_raw_aligned[:, 1], z_raw_aligned[:, 2]))])
         total_gradient = np.vstack([total_gradient, np.column_stack((gradient_xy[:, 0], gradient_xy[:, 1], gradient_xy[:, 2], gradient_xy[:, 3]))])
     
     if 0:
@@ -185,7 +282,7 @@ def main(config_path, img_folder):
         ax.set_xlabel('X (mm)')
         ax.set_ylabel('Y (mm)')
         ax.set_zlabel('Z (mm)')
-        ax.set_title('3D Point Cloud Visualization gradient')
+        ax.set_title('z_raw 3D Point Cloud Visualization gradient')
         plt.show()
 
         # fig = plt.figure()
@@ -215,13 +312,13 @@ def main(config_path, img_folder):
     total_cloud_point[:,0] = np.round(total_cloud_point[:,0])
     total_cloud_point[:,1] = np.round(total_cloud_point[:,1])
     #fitted_points = post_process(total_cloud_point, debug=0)
-    fitted_points = post_process_with_grad(img_folder, 0)
+    fitted_points = post_process_with_grad(img_folder, n_cam, 1)
     #fitted_points = total_cloud_point
     #align_fitted, _ = align2ref(fitted_points)
     align_fitted = fitted_points
     
     write_point_cloud(os.path.join(img_folder, 'cloudpoint.txt'), np.round(align_fitted[:, 0]-np.mean(align_fitted[:, 0])), np.round(align_fitted[:, 1]-np.mean(align_fitted[:, 1])), align_fitted[:,2]-np.min(align_fitted[:,2]))
-    if 0:
+    if 1:
         fig = plt.figure()
         ax = fig.add_subplot(111, projection='3d')
 
@@ -230,6 +327,10 @@ def main(config_path, img_folder):
         y_vals = np.round(fitted_points[:, 1]-np.mean(fitted_points[:, 1]))
         z_vals = align_fitted[:,2]-np.min(align_fitted[:,2])
 
+        x_vals = fitted_points[:, 0]
+        y_vals = fitted_points[:, 1]
+        z_vals = fitted_points[:, 2]
+
         # 绘制3D点云
         ax.scatter(x_vals, y_vals, z_vals, c=z_vals, cmap='viridis', marker='o')
 
@@ -237,7 +338,7 @@ def main(config_path, img_folder):
         ax.set_xlabel('X (mm)')
         ax.set_ylabel('Y (mm)')
         ax.set_zlabel('Z (mm)')
-        ax.set_title('3D Point Cloud Visualization')
+        ax.set_title('post 3D Point Cloud Visualization')
         plt.show()
     
     post_process_end = time.time()
@@ -249,7 +350,7 @@ def main(config_path, img_folder):
     point_cloud_path = os.path.join(img_folder, 'cloudpoint.txt')
     json_path = os.path.join(img_folder, 'result.json')
     theta_notch = 0
-    get_eval_result(point_cloud_path, json_path, theta_notch, 0)
+    #get_eval_result(point_cloud_path, json_path, theta_notch, 0)
     eval_end = time.time()
     print(f"   完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
     print(f"   耗时: {eval_end - eval_start:.2f} 秒")
@@ -257,14 +358,14 @@ def main(config_path, img_folder):
     return True
 
 
-
-
-
 if __name__ == '__main__':
-    config_path = 'config\\cfg_3freq_wafer.json'
+    config_path = 'config\\cfg_3freq_wafer_matlab.json'
     #img_folder = 'D:\\data\\four_cam\\1008_storage\\pingjing_20241017092315228\\'   #'D:\\data\\four_cam\\betone_1011\\20241011142348762-1'
    # img_folder = 'D:\\huchao\\inspect_server_202409241013_py\\storage\\20241023193453708\\'
-    img_folder = 'D:\\data\\four_cam\\betone_1011\\20241011144821292-4\\'
+    img_folder = 'D:\\data\\one_cam\\pingjing_20241024154405250\\'
+    #img_folder = 'D:\\data\\one_cam\\betone-20241025095352783\\'
+    #img_folder = 'D:\\data\\one_cam\\20241025113857041-neg\\'
+    #'
     json_path = os.path.join(img_folder, 'result.json')
 
     # 20241011142348762-1  20241011142901251-2   20241011143925746-3   20241011144821292-4 
@@ -281,7 +382,5 @@ if __name__ == '__main__':
 
     pmdstart(config_path, img_folder)
     #fitted_points = post_process_with_grad(img_folder, 1)
-
-    
     
     

src/calibration/calibrate.py

@@ -7,6 +7,8 @@ def map_screen_to_world(data_screen, date_world):
     A = date_world['camera_matrix']
     Rw2c = date_world['rotation_matrix']
     Tw2c = date_world['translation_vector']
+    # print('Rv2c = ', Rv2c)
+    # print('Tv2c = ', Tv2c)
 
     # print(Rv2c.shape, Tv2c.shape, A.shape, Rw2c.shape, Tw2c.shape)
     # Rv2c = np.array([[0.965149963314468, -0.0425788059963296, 0.258210366937518],

src/calibration/get_camera_params.py

@@ -169,6 +169,9 @@ def calibrate_world(calibration_images_path, cam_id, world_chessboard_size=(6,6)
     mean_error = total_error / len(all_object_points)
     #print(f"Mean re-projection error: {mean_error}")
 
+    print('rvecs[0] = ', rvecs[0])
+    #print('rvecs[0] = ', (rvecs[0]))
+
     # 获得第一张图片的旋转矩阵和平移向量
     R = cv2.Rodrigues(rvecs[0])[0]  # 转换为旋转矩阵
     T = tvecs[0]

src/recons.py

@@ -311,7 +311,7 @@ def poisson_recons(x, y, gx, gy, debug=False):
 
 
 
-def reconstruction_cumsum(world_points, camera_points, screen_points, save_path, cam_id, debug=False, smooth=False, align=True, denoise=True):
+def reconstruction_cumsum(world_points, camera_points, screen_points, debug=False, smooth=False, align=True, denoise=True):
 
     x = world_points[:, 0]
     y = world_points[:, 1]
@@ -338,7 +338,6 @@ def reconstruction_cumsum(world_points, camera_points, screen_points, save_path,
 
     #print('normal_vector = ', normal_vector[0])
 
-
     # 计算法向量在X方向的梯度
     gradient_x = normal_vector[0]
 
@@ -347,6 +346,9 @@ def reconstruction_cumsum(world_points, camera_points, screen_points, save_path,
     #gradient_x = np.zeros(gradient_x.shape)
     #gradient_y = np.zeros(gradient_x.shape)
 
+    print('abs mean (gradient_x) = ', np.mean(abs(gradient_x)))
+    print('abs mean (gradient_y) = ', np.mean(abs(gradient_y)))
+
     #gradient_x = gradient_x - np.mean(gradient_x)
     #gradient_y = gradient_y - np.mean(gradient_y)
     # import pdb
@@ -355,18 +357,13 @@ def reconstruction_cumsum(world_points, camera_points, screen_points, save_path,
     z = np.zeros(len(x))
     z[0] = 0  # 初始点的z
     
-    #z_raw = poisson_recons(x, y, gradient_x, gradient_y)
-    z_raw = poisson_recons_with_smoothed_gradients(x, y, gradient_x, gradient_y)
+    z_raw = poisson_recons(x, y, gradient_x, gradient_y)
+    #z_raw = poisson_recons_with_smoothed_gradients(x, y, gradient_x, gradient_y)
+    #z_raw = z
 
     #print('gradient_x  = ', x.shape, y.shape, gradient_x.shape, gradient_y.shape)
-
     gradient_xy = np.column_stack((x, y, gradient_x, gradient_y))
     
-       # z_raw = least_squares_recons(x, y, gradient_x, gradient_y)
-        # import pdb; pdb.set_trace()
-    #print('-- point cloud: mean - {}; std - {}'.format(np.nanmean(z), np.nanstd(z)))
-    
-        #print('-- nanmean after smoothing: mean - {}; std - {}'.format(np.nanmean(smoothed_points[:, 2]), np.nanstd(smoothed_points[:, 2])))
     if align:
         #points_aligned, rotation_matrix = align2ref(smoothed_points)
         points_aligned, rotation_matrix = align2ref(np.column_stack((x, y, z_raw)))
@@ -389,31 +386,6 @@ def reconstruction_cumsum(world_points, camera_points, screen_points, save_path,
         # points_denoise = remove_noise(x, y, z, a, b, c, thresh=0.01)
 
     if debug:
-        #fig, axes = plt.subplots(1, 2, figsize=(12, 6))
-
-        # # 第一个子图
-        # cax0 = axes[0].scatter(x, y, c=gradient_x, cmap='viridis')
-        # axes[0].set_title('Visualization of x gradient')
-        # axes[0].set_xlabel('X Axis')
-        # axes[0].set_ylabel('Y Axis')
-        # fig.colorbar(cax0, ax=axes[0])
-
-        # cax1 = axes[1].scatter(x, y, c=gradient_y, cmap='viridis')
-        # axes[1].set_title('Visualization of y gradient')
-        # axes[1].set_xlabel('X Axis')
-        # axes[1].set_ylabel('Y Axis')
-        # fig.colorbar(cax1, ax=axes[1])
-        # #plt.savefig(os.path.join(save_path, "gradient.png"))
-        # # 调整子图之间的间距
-        # plt.tight_layout()
-        # # 显示图像
-        # plt.show()
-        # fig = plt.figure()
-        # ax = fig.add_subplot(111, projection='3d')
-        # # 绘制3D点云
-        # ax.scatter(x, y, z, c=z, cmap='viridis', marker='o')
-        # plt.show()
-
         fig = plt.figure()
         ax = fig.add_subplot(111, projection='3d')
 
@@ -431,27 +403,9 @@ def reconstruction_cumsum(world_points, camera_points, screen_points, save_path,
         ax.set_zlabel('Z (mm)')
         ax.set_title('raw 3D Point Cloud Visualization')
 
-        fig = plt.figure()
-        ax = fig.add_subplot(111, projection='3d')
-
-        # 提取 x, y, z 坐标
-        x_vals = x
-        y_vals = y
-        z_vals = points_smoothed[:,2]
-
-        # 绘制3D点云
-        ax.scatter(x_vals, y_vals, z_vals, c=z_vals, cmap='viridis', marker='o', s=1)
-
-        # 设置轴标签和标题
-        ax.set_xlabel('X (mm)')
-        ax.set_ylabel('Y (mm)')
-        ax.set_zlabel('Z (mm)')
-        ax.set_title('smooth 3D Point Cloud Visualization')
         plt.show()
-
-   
-    
-    return np.column_stack((x, y, z_raw)), points_aligned, points_smoothed, points_denoise, gradient_xy
+ 
+    return np.column_stack((x, y, z_raw)), gradient_xy
 
 
 

src/utils.py

@@ -369,6 +369,7 @@ def get_world_points(contour_points, world_mat_contents, cam_id, grid_spacing, d
     Tc2w = -np.dot(np.linalg.inv(Rw2c), Tw2c)
     world_points_contours = affine_img2world(contour_points, A, Rc2w, Tc2w, d)
     
+    
     world_points_fit,  world_points_boundary = fit_sector(world_points_contours, grid_spacing, erosion_pixels)
     _,  world_points_boundary_3 = fit_sector(world_points_contours, grid_spacing, 3)
     # assert np.abs(world_contour_points[:, 0].max()-world_contour_points[:, 0].min()-(world_contour_points[:, 1].max()-world_contour_points[:, 1].min())) < 1
@@ -496,10 +497,10 @@ def get_screen_points(point_data, x_phase_unwrapped, y_phase_unwrapped, screen_p
                             0])
     abs_phasepoint_world = np.dot(Rs2w, arrow_abs_mm) + Ts2w.T
     
-    # print('abs_phasepoint 绝对相位值:')
-    # print(abs_phasepoint)
-    # print('abs_phasepoint_world 绝对相位值对应的世界坐标系的位置:')
-    # print(abs_phasepoint_world)
+    print('abs_phasepoint 绝对相位值:')
+    print(abs_phasepoint)
+    print('abs_phasepoint_world 绝对相位值对应的世界坐标系的位置:')
+    print(abs_phasepoint_world)
 
     # print('----------------------------------------')
 
@@ -582,7 +583,10 @@ def get_white_mask(white_path, bin_thresh=15, debug=False):
         max_index, max_contour = max(enumerate(valid_contours), key=lambda x: cv2.boundingRect(x[1])[2] * cv2.boundingRect(x[1])[3]) 
         cv2.drawContours(mask, [max_contour], -1, (255, 255, 255), thickness=cv2.FILLED)
     if debug:
-       
+        
+        cv2.namedWindow('gray',0)
+        cv2.imshow('gray', gray)
+        cv2.waitKey(0)
    
         cv2.namedWindow('mask',0)
         cv2.imshow('mask', mask)
@@ -801,12 +805,12 @@ def curve_fit_gaussian(points):
     # 提取拟合参数
     A, x0, y0, sigma = popt
     
-    #print(f"拟合参数: A={A}, x0={x0}, y0={y0}, sigma={sigma}")
-    x_mean = np.mean(x)
-    y_mean = np.mean(y)
-    x0 = x_mean
-    y0 = y_mean
-    sigma = 2 * sigma
+    print(f"拟合参数: A={A}, x0={x0}, y0={y0}, sigma={sigma}")
+    # x_mean = np.mean(x)
+    # y_mean = np.mean(y)
+    # x0 = x_mean
+    # y0 = y_mean
+    # sigma = 2 * sigma
     #print(f"拟合参数1: A={A}, x0={x0}, y0={y0}, sigma={sigma}")
     
 
@@ -1240,9 +1244,9 @@ def poisson_reconstruct(x_grad_expand, y_grad_expand):
 
 
 
-def post_process_with_grad(img_folder, debug=False):
+def post_process_with_grad(img_folder, n_cam, debug=False):
     total_point_grad = np.empty((0, 4))
-    for i in range(4):
+    for i in range(n_cam):
         txt_file = img_folder + str(i) + '_gradient.txt'
         total_point_grad = np.vstack([total_point_grad, np.loadtxt(os.path.join(txt_file), delimiter=',')])
      
@@ -1255,27 +1259,45 @@ def post_process_with_grad(img_folder, debug=False):
     x = total_point_grad[:,0]
     y = total_point_grad[:,1]
 
-    x_gradient = total_point_grad[:,2]
-    y_gradient = total_point_grad[:,3]
+    x_gradient = total_point_grad[:,2] 
+    y_gradient = total_point_grad[:,3] 
+
+    mean_x_grad = np.mean((x_gradient))
+    mean_y_grad = np.mean((y_gradient))
 
-    mean_x_grad = np.mean(x_gradient)
-    mean_y_grad = np.mean(y_gradient)
+    sum_x = np.sum(abs(x_gradient))
+    sum_y = np.sum(abs(y_gradient))
+
+    print('sum x grad = ', sum_x)
+    print('sum y grad = ', sum_y)
+
+    print('mean x grad = ', mean_x_grad)
+    print('mean y grad = ', mean_y_grad)
 
     x_gradient =  (x_gradient - mean_x_grad)
     y_gradient =  (y_gradient - mean_y_grad)
     
+    print('max min')
+    print(max(x_gradient), min(x_gradient))
+    print(max(y_gradient), min(y_gradient))
+    print('max min')
+    
     # 使用最小外接矩形，并填充梯度
     x_grad_expand, y_grad_expand, x_expand, y_expand = expand_to_rectangle_with_unit_spacing(x, y, x_gradient, y_gradient)   
     #print(x_grad_expand.shape, y_grad_expand.shape, x_expand.shape, y_expand.shape)
     # 调用 Southwell 重建函数
     #Z_reconstructed = southwell_integrate_smooth(x_grad_expand, y_grad_expand, x_expand, y_expand, (0, 0))
 
-    Z_reconstructed = poisson_reconstruct(x_grad_expand, y_grad_expand)
+    Z_reconstructed = np.cumsum(x_grad_expand, 1)*1 + np.cumsum(y_grad_expand, 0)*1
+    print('min_z =', np.min(Z_reconstructed))
+
+    #Z_reconstructed = poisson_reconstruct(x_grad_expand, y_grad_expand)
     
     z1_lookup = {(x1, y1): z1 for x1, y1, z1 in np.column_stack((x_expand.reshape(-1,1), y_expand.reshape(-1,1), Z_reconstructed.reshape(-1,1)))}
 
     z1_values = np.array([z1_lookup.get((x, y), np.nan) for x, y in np.column_stack((x, y))])
-    fitted_points = post_process(np.column_stack((x,y,z1_values)), debug=0)
+    #fitted_points = post_process(np.column_stack((x,y,z1_values)), debug=0)
+
 
     # 可视化梯度场
     if debug:
@@ -1283,12 +1305,12 @@ def post_process_with_grad(img_folder, debug=False):
         ax = fig.add_subplot(111, projection='3d')
 
         #提取 x, y, z 坐标
-        x_vals = fitted_points[:,0]
-        y_vals = fitted_points[:,1]
-        z_vals = fitted_points[:,2]
-        # x_vals = x
-        # y_vals = y
-        # z_vals = z1_values
+        # x_vals = fitted_points[:,0]
+        # y_vals = fitted_points[:,1]
+        # z_vals = fitted_points[:,2]
+        x_vals = x_expand
+        y_vals = y_expand
+        z_vals = x_grad_expand
         # 绘制3D点云
         ax.scatter(x_vals, y_vals, z_vals, c=z_vals, cmap='viridis', marker='o')
 
@@ -1299,7 +1321,8 @@ def post_process_with_grad(img_folder, debug=False):
         ax.set_title('post_process 3D Point Cloud Visualization')
         plt.show()
     
-    return fitted_points
+    #return fitted_points
+    return np.column_stack((x,y,z1_values))
 
 
 
@@ -1420,7 +1443,7 @@ def find_best_circle(image, min_thresh, max_thresh):
 
 
 
-def find_notch(white_path, cam_num):
+def find_notch(white_path, cam_num, debug=0):
 
     image = cv2.imread(white_path)
     if cam_num == 1:
@@ -1450,15 +1473,16 @@ def find_notch(white_path, cam_num):
         angle_rad = np.arctan2(min_pt[1] - center_y, min_pt[0] - center_x)
 
         thresh_img = cv2.drawContours(background, [max_contour_ellipse], -1, (255, 255, 255), thickness=cv2.FILLED)
-    
-        cv2.circle(image, min_pt, 50, (0,255,255), 2)
-        #print(min_pt)
-        cv2.namedWindow('thresh_img', 0)
-        cv2.imshow('thresh_img', thresh_img)
-        cv2.namedWindow('contours', 0)
-        cv2.imshow('contours', image)
-    
-        cv2.waitKey(0)
+
+        if debug:
+            cv2.circle(image, min_pt, 50, (0,255,255), 2)
+            #print(min_pt)
+            cv2.namedWindow('thresh_img', 0)
+            cv2.imshow('thresh_img', thresh_img)
+            cv2.namedWindow('contours', 0)
+            cv2.imshow('contours', image)
+        
+            cv2.waitKey(0)
         return angle_rad, thresh_img
     else:
         return 0
@@ -1870,7 +1894,15 @@ def img_diff():
         cv2.waitKey(0)
         
         
-
+def mat2pkl():
+    from scipy.io import loadmat, savemat
+    mat_path = 'D:\\code\\pmd-python\\config_1\\screen_params.mat'
+    pkl_path = 'D:\\code\\pmd-python\\config\\screen_params_1.pkl'
+    matdata = loadmat(mat_path)
+    print(matdata)
+    print([matdata])
+    with open(pkl_path, 'wb') as pkl_file:
+        pickle.dump([matdata], pkl_file)