Approximate a 2D surface.
(
data: Array,
u: Array,
v: Array,
uorder: int,
vorder: int,
uknots: int | Array = 50,
vknots: int | Array = 50,
uperiodic: bool = False,
vperiodic: bool = False,
penalty: int = 3,
lam: float = 0,
)
| 1811 | |
| 1812 | |
| 1813 | def approximate2D( |
| 1814 | data: Array, |
| 1815 | u: Array, |
| 1816 | v: Array, |
| 1817 | uorder: int, |
| 1818 | vorder: int, |
| 1819 | uknots: int | Array = 50, |
| 1820 | vknots: int | Array = 50, |
| 1821 | uperiodic: bool = False, |
| 1822 | vperiodic: bool = False, |
| 1823 | penalty: int = 3, |
| 1824 | lam: float = 0, |
| 1825 | ) -> Surface: |
| 1826 | """ |
| 1827 | Approximate a 2D surface. |
| 1828 | """ |
| 1829 | |
| 1830 | # process the knots |
| 1831 | uknots_ = uknots if isinstance(uknots, Array) else np.linspace(0, 1, uknots) |
| 1832 | vknots_ = vknots if isinstance(vknots, Array) else np.linspace(0, 1, vknots) |
| 1833 | |
| 1834 | # create the design matrix |
| 1835 | C = designMatrix2D( |
| 1836 | u, v, uorder, vorder, uknots_, vknots_, uperiodic, vperiodic |
| 1837 | ).csc() |
| 1838 | |
| 1839 | # handle penalties if requested |
| 1840 | if lam: |
| 1841 | # construct the penalty grid |
| 1842 | up, vp = uniformGrid(uknots_, vknots_, uorder, vorder, uperiodic, vperiodic) |
| 1843 | |
| 1844 | # construct the derivative matrices |
| 1845 | penalties = penaltyMatrix2D( |
| 1846 | up, vp, uorder, vorder, penalty, uknots_, vknots_, uperiodic, vperiodic, |
| 1847 | ) |
| 1848 | |
| 1849 | # augment the design matrix |
| 1850 | tmp = [comb(penalty, i) * penalties[i].csc() for i in range(penalty + 1)] |
| 1851 | Lu = uknots_[-1] - uknots_[0] # v length of the parametric domain |
| 1852 | Lv = vknots_[-1] - vknots_[0] # u length of the parametric domain |
| 1853 | P = Lu * Lv / len(up) * sp.vstack(tmp) |
| 1854 | |
| 1855 | CtC = C.T @ C + lam * P.T @ P |
| 1856 | else: |
| 1857 | CtC = C.T @ C |
| 1858 | |
| 1859 | # solve normal equations |
| 1860 | CtC = DM(CtC) |
| 1861 | ldl = Linsol("mumps", "mumps", CtC.sparsity()) |
| 1862 | pts = ldl.solve(CtC, C.T @ data).toarray() |
| 1863 | |
| 1864 | # construct the result |
| 1865 | rv = Surface( |
| 1866 | pts.reshape((len(uknots_) - int(uperiodic), len(vknots_) - int(vperiodic), 3)), |
| 1867 | uknots_, |
| 1868 | vknots_, |
| 1869 | uorder, |
| 1870 | vorder, |