|
21 | 21 | from .vec_utils import unit_vector, vec_norm, vec_set |
22 | 22 |
|
23 | 23 |
|
| 24 | +def is_singular(matrix): |
| 25 | + rank = np.linalg.matrix_rank(matrix) |
| 26 | + return rank < matrix.shape[0] |
| 27 | + |
24 | 28 | @njit |
25 | 29 | def skew_midpoint( |
26 | 30 | vert1: np.ndarray, direct1: np.ndarray, vert2: np.ndarray, direct2: np.ndarray |
@@ -578,13 +582,29 @@ def orient( |
578 | 582 | # if np.any(np.isnan(X)): |
579 | 583 | # pdb.set_trace() |
580 | 584 |
|
581 | | - XPX = np.linalg.inv(np.dot(X[:, :numbers].T, X[:, :numbers])) |
| 585 | + XTX = np.dot(X[:, :numbers].T, X[:, :numbers]) |
| 586 | + if is_singular(XTX): |
| 587 | + XPX = np.linalg.pinv(XTX) |
| 588 | + else: |
| 589 | + XPX = np.linalg.inv(XTX) |
| 590 | + |
| 591 | + |
| 592 | + # XPX = np.linalg.inv() |
| 593 | + |
| 594 | + # def invert_singular_matrix(m): |
| 595 | + # a, b = m.shape |
| 596 | + # if a != b: |
| 597 | + # raise ValueError("Only square matrices are invertible.") |
| 598 | + # identity_matrix = np.eye(a, a) |
| 599 | + # return np.linalg.lstsq(m, identity_matrix)[0] |
582 | 600 |
|
583 | 601 |
|
584 | 602 | # import pdb; pdb.set_trace() |
585 | | - for i in range(numbers): |
586 | | - # print(f"{i=}, {np.sqrt(XPX[i][i]) = }") |
587 | | - sigmabeta[i] = sigmabeta[NPAR] * np.sqrt(XPX[i][i]) |
| 603 | + # for i in range(numbers): |
| 604 | + # # print(f"{i=}, {np.sqrt(XPX[i][i]) = }") |
| 605 | + # sigmabeta[i] = sigmabeta[NPAR] * np.sqrt(XPX[i][i]) |
| 606 | + |
| 607 | + sigmabeta[:numbers] = sigmabeta[NPAR]*np.sqrt(np.diag(XPX)) |
588 | 608 |
|
589 | 609 | if stopflag: |
590 | 610 | cal.update_rotation_matrix() |
@@ -843,13 +863,16 @@ def full_calibration( |
843 | 863 | """ |
844 | 864 | err_est = np.empty((NPAR + 1), dtype=np.float64) |
845 | 865 |
|
846 | | - # convert numpy array to list of Target objects |
847 | | - targs = [Target() for _ in img_pts] |
| 866 | + if isinstance(img_pts, np.ndarray): |
| 867 | + # convert numpy array to list of Target objects |
| 868 | + targs = [Target() for _ in img_pts] |
848 | 869 |
|
849 | | - for ptx, pt in enumerate(img_pts): |
850 | | - targs[ptx].x = pt[0] |
851 | | - targs[ptx].y = pt[1] |
852 | | - targs[ptx].pnr = ptx |
| 870 | + for ptx, pt in enumerate(img_pts): |
| 871 | + targs[ptx].x = pt[0] |
| 872 | + targs[ptx].y = pt[1] |
| 873 | + targs[ptx].pnr = ptx |
| 874 | + else: |
| 875 | + targs = img_pts |
853 | 876 |
|
854 | 877 | residuals = orient(cal, cparam, len(ref_pts), ref_pts, targs, orient_par, err_est, dm=dm, drad=drad) |
855 | 878 |
|
|
0 commit comments