Multiplane calibration with additional parameters (#71)

* added scipy.optimize.minimize option. fails on plot

* added editing addpar files

* fixed saving addpar and ori files

* working with k1 only, need to make it generic

* fully working version for multiplane calibration with provided values and additional parameters

milestone achieved for the proPTV, MyPTV, OpenPTV benchmarking
- use markers as xy, XYZ and convert them to .crd and .fix files
- use multiplane combine to read those files in combined mode
- use fast C calibration for external parameters, back focal point (cc, xh, yh)
- use scipy.optimize.minimize to the radial distortion, decentering (p1,p2) and affine transformation (scale, shear)
\

pyptv/calibration_gui.py

@@ -32,7 +32,7 @@
 
 # from chaco.tools.simple_zoom import SimpleZoom
 from pyptv.text_box_overlay import TextBoxOverlay
-from pyptv.code_editor import codeEditor
+from pyptv.code_editor import oriEditor, addparEditor
 from pyptv.quiverplot import QuiverPlot
 
 
@@ -47,6 +47,7 @@
 
 from pyptv import ptv, parameter_gui, parameters as par
 
+from scipy.optimize import minimize
 
 # -------------------------------------------
 class ClickerTool(ImageInspectorTool):
@@ -340,6 +341,7 @@ class CalibrationGUI(HasTraits):
     button_checkpoint = Button()
     button_ap_figures = Button()
     button_edit_ori_files = Button()
+    button_edit_addpar_files = Button()
     button_test = Button()
 
     # ---------------------------------------------------
@@ -480,6 +482,11 @@ def __init__(self, active_path: Path):
                         label="Edit ori files",
                         show_label=False,
                     ),
+                    Item(
+                        name="button_edit_addpar_files",
+                        label="Edit addpar files",
+                        show_label=False,
+                    ),
                     show_left=False,
                 ),
             ),
@@ -872,8 +879,9 @@ def _button_fine_orient_fired(self):
             op.shear,
         ]
 
+        # set flags for cc, xh, yh only
         flags = []
-        for name, op_name in zip(names, op_names):
+        for name, op_name in zip(names[:3], op_names[:3]):
             if op_name == 1:
                 flags.append(name)
 
@@ -968,24 +976,99 @@ def _button_fine_orient_fired(self):
             else:
                 targs = self.sorted_targs[i_cam]
 
-            try:
-                residuals, targ_ix, err_est = full_calibration(
-                    self.cals[i_cam],
-                    self.cal_points["pos"],
-                    targs,
-                    self.cpar,
-                    flags,
-                )
-            except BaseException as exc:
-                raise ValueError("full calibration failed\n") from exc
-            # save the results
+            # try:
+            print(f"First calibrate only external and flags: {flags} \n")
+            residuals, targ_ix, err_est = full_calibration(
+                self.cals[i_cam],
+                self.cal_points["pos"],
+                targs,
+                self.cpar,
+                flags,
+            )
+
+            # this chunk optimizes for radial distortion
+            if np.any(op_names[3:6]):
+                sol = minimize(self._residuals_k,
+                               self.cals[i_cam].get_radial_distortion(), 
+                               args=(self.cals[i_cam], 
+                                     all_known, 
+                                     all_detected, 
+                                     self.cpar
+                                     ), 
+                                     method='Nelder-Mead', 
+                                     tol=1e-11,
+                                     options={'disp':True},
+                                     )
+                radial = sol.x 
+                self.cals[i_cam].set_radial_distortion(radial)
+            else:
+                radial = self.cals[i_cam].get_radial_distortion()
+
+            if np.any(op_names[5:8]):
+                # now decentering
+                sol = minimize(self._residuals_p,
+                               self.cals[i_cam].get_decentering(), 
+                               args=(self.cals[i_cam], 
+                                     all_known, 
+                                     all_detected, 
+                                     self.cpar
+                                     ), 
+                                     method='Nelder-Mead', 
+                                     tol=1e-11,
+                                     options={'disp':True},
+                                     )
+                decentering = sol.x 
+                self.cals[i_cam].set_decentering(decentering)
+            else:
+                decentering = self.cals[i_cam].get_decentering()
+
+            if np.any(op_names[8:]):
+                # now affine
+                sol = minimize(self._residuals_s,
+                               self.cals[i_cam].get_affine(), 
+                               args=(self.cals[i_cam], 
+                                     all_known, 
+                                     all_detected, 
+                                     self.cpar
+                                     ), 
+                                     method='Nelder-Mead', 
+                                     tol=1e-11,
+                                     options={'disp':True},
+                                     )
+                affine = sol.x 
+                self.cals[i_cam].set_affine_trans(affine)
+
+            else:
+                affine = self.cals[i_cam].get_affine()
+
+
+            # Now project and estimate full residuals
+            self._project_cal_points(i_cam)
+
+            residuals = self._residuals_combined(
+                            np.r_[radial, decentering, affine],
+                            self.cals[i_cam], 
+                            all_known, 
+                            all_detected, 
+                            self.cpar
+                            )
+
+            residuals /= 100
+
+            targ_ix = np.arange(len(all_detected))
+
+            # save the results from self.cals[i_cam]
             self._write_ori(i_cam, addpar_flag=True)
 
-            # Plot the output
-            # self.reset_plots()
+            # x, y = [], []
+            # for r, t in zip(residuals, targ_ix):
+            #     if t != -999:
+            #         pos = targs[t].pos()
+            #         x.append(pos[0])
+            #         y.append(pos[1])
 
             x, y = [], []
-            for r, t in zip(residuals, targ_ix):
+            for t in targ_ix:
                 if t != -999:
                     pos = targs[t].pos()
                     x.append(pos[0])
@@ -1012,6 +1095,78 @@ def _button_fine_orient_fired(self):
 
         self.status_text = "Orientation finished."
 
+    # def _error_function(self, x, cal, XYZ, xy, cpar):
+    #     """Error function for scipy.optimize.minimize.
+
+    #     Args:
+    #         x (array-like): Array of parameters.
+    #         cal (Calibration): Calibration object.
+    #         XYZ (array-like): 3D coordinates.
+    #         xy (array-like): 2D image coordinates.
+    #         cpar (CPar): Camera parameters.
+
+    #     Returns:
+    #         float: Error value.
+    #     """
+    #     residuals = self._residuals_radial(x, cal, XYZ, xy, cpar)
+    #     return np.sum(residuals**2)
+
+    def _residuals_k(self, x, cal, XYZ, xy, cpar):
+        """Residuals due to radial distortion
+
+        Args:
+            x (array-like): Array of parameters.
+            cal (Calibration): Calibration object.
+            XYZ (array-like): 3D coordinates.
+            xy (array-like): 2D image coordinates.
+            cpar (CPar): Camera parameters.
+
+        Returns:
+            residuals: Distortion in pixels
+        """
+
+        cal.set_radial_distortion(x)
+        targets = convert_arr_metric_to_pixel(
+            image_coordinates(XYZ, cal, cpar.get_multimedia_params()),
+            cpar
+        )
+        residuals = xy[:,1:] - targets
+        return np.sum(residuals**2)
+
+    def _residuals_p(self, x, cal, XYZ, xy, cpar):
+        """Residuals due to decentering """
+        cal.set_decentering(x)
+        targets = convert_arr_metric_to_pixel(
+            image_coordinates(XYZ, cal, cpar.get_multimedia_params()),
+            cpar
+        )
+        residuals = xy[:,1:] - targets
+        return np.sum(residuals**2)
+
+    def _residuals_s(self, x, cal, XYZ, xy, cpar):
+        """Residuals due to decentering """
+        cal.set_affine_trans(x)
+        targets = convert_arr_metric_to_pixel(
+            image_coordinates(XYZ, cal, cpar.get_multimedia_params()),
+            cpar
+        )
+        residuals = xy[:,1:] - targets
+        return np.sum(residuals**2)    
+
+    def _residuals_combined(self, x, cal, XYZ, xy, cpar):
+        """Combined residuals  """
+
+        cal.set_radial_distortion(x[:3])
+        cal.set_decentering(x[3:5])
+        cal.set_affine_trans(x[5:])
+
+        targets = convert_arr_metric_to_pixel(
+            image_coordinates(XYZ, cal, cpar.get_multimedia_params()),
+            cpar
+        )
+        residuals = xy[:,1:] - targets
+        return residuals           
+
     def _write_ori(self, i_cam, addpar_flag=False):
         """Writes ORI and ADDPAR files for a single calibration result
         of i_cam
@@ -1116,9 +1271,14 @@ def reset_show_images(self):
             cam._plot.request_redraw()
 
     def _button_edit_ori_files_fired(self):
-        editor = codeEditor(path=self.par_path)
+        editor = oriEditor(path=self.par_path)
+        editor.edit_traits(kind="livemodal")
+
+    def _button_edit_addpar_files_fired(self):
+        editor = addparEditor(path=self.par_path)
         editor.edit_traits(kind="livemodal")
 
+
     def drawcross(self, str_x, str_y, x, y, color1, size1, i_cam=None):
         """
 
@@ -1174,6 +1334,7 @@ def _read_cal_points(self):
 
         with open(self.calParams.fixp_name, 'r') as file:
             first_line = file.readline()
+            print(first_line)
             if ',' in first_line:
                 delimiter=','
             elif '\t' in first_line:
@@ -1198,6 +1359,7 @@ def _read_cal_points(self):
 
     if len(sys.argv) == 1:
         active_path = Path("../test_cavity/parametersRun1")
+        active_path = Path("/home/user/Downloads/rbc300/parametersMultiPlane")
     else:
         active_path = Path(sys.argv[0])
 

pyptv/code_editor.py

@@ -28,41 +28,47 @@ def get_path(filename):
     )
 
 
-def get_code(path):
-    f = open(path, "r")
+def get_code(path: Path):
+    """ Read the code from the file """
+
+    # print(f"Read from {path}: {path.exists()}")
+    with open(path, "r", encoding="utf-8") as f:    
+        retCode = f.read()
+
+    # print(retCode)
 
-    retCode = f.read()
-    f.close()
     return retCode
 
 
-class oriEditor(HasTraits):
-    file_Path = File
-    ori_Code = Code()
-    ori_Save = Button(label="Save")
+class codeEditor(HasTraits):
+    file_Path = Path
+    _Code = Code()
+    save_button = Button(label="Save")
     buttons_group = Group(
-        Item(name="file_Path", style="simple", show_label=False, width=0.3),
-        Item(name="ori_Save", show_label=False),
+        Item(name="file_Path", style="simple", show_label=True, width=0.3),
+        Item(name="save_button", show_label=True),
         orientation="horizontal",
     )
     traits_view = View(
         Group(
-            Item(name="ori_Code", show_label=False, height=300, width=650),
+            Item(name="_Code", show_label=False, height=300, width=650),
             buttons_group,
         )
     )
 
-    def _ori_Save_fired(self, filename, code):
-        f = open(self.file_Path, "w")
-        f.write(self.ori_Code)
-        f.close()
-
-    def __init__(self, file_path):
+    def _save_button_fired(self):
+        with open(self.file_Path, "w", encoding="utf-8") as f:
+            # print(f"Saving to {self.file_Path}")
+            # print(f"Code: {self._Code}")
+            f.write(self._Code)
+
+        print(f"Saved to {self.file_Path}")
+
+    def __init__(self, file_path: Path):
         self.file_Path = file_path
-        self.ori_Code = get_code(file_path)
-
+        self._Code = get_code(file_path)
 
-class codeEditor(HasTraits):
+class oriEditor(HasTraits):
 
     # number of images
     n_img = Int()
@@ -86,6 +92,47 @@ class codeEditor(HasTraits):
         title="Camera's orientation files",
     )
 
+    def __init__(self, path: Path):
+        """ Initialize by reading parameters and filling the editor windows """
+        # load ptv_par
+        ptvParams = par.PtvParams(path=path)
+        ptvParams.read()
+        self.n_img = ptvParams.n_img
+
+        # load cal_ori
+        calOriParams = par.CalOriParams(self.n_img)
+        calOriParams.read()
+
+        for i in range(self.n_img):
+            self.oriEditors.append(
+                codeEditor(Path(calOriParams.img_ori[i]))
+            )
+
+
+class addparEditor(HasTraits):
+
+    # number of images
+    n_img = Int()
+
+    addparEditors = List
+
+    # view
+    traits_view = View(
+        Item(
+            "addparEditors",
+            style="custom",
+            editor=ListEditor(
+                use_notebook=True,
+                deletable=False,
+                dock_style="tab",
+                page_name=".file_Path",
+            ),
+            show_label=False,
+        ),
+        buttons=["Cancel"],
+        title="Camera's additional parameters files",
+    )
+
     def __init__(self, path):
         """ Initialize by reading parameters and filling the editor windows """
         # load ptv_par
@@ -98,6 +145,6 @@ def __init__(self, path):
         calOriParams.read()
 
         for i in range(self.n_img):
-            self.oriEditors.append(
-                oriEditor(calOriParams.img_ori[i])
+            self.addparEditors.append(
+                codeEditor(Path(calOriParams.img_ori[i].replace('ori', 'addpar')))
             )