diff --git a/morphops/io.py b/morphops/io.py index e6808e6..ca05050 100644 --- a/morphops/io.py +++ b/morphops/io.py @@ -120,7 +120,7 @@ def write_dta(filename, lmk_sets, names=[]): # Write the names. Missing names are populated as 'InsertName{ID}', # where ID goes from [len(names) + 1, n + 1). rem_names_ids = np.arange(len(names)+1,n+1).astype(str) - rem_names = np.core.char.add('InsertName', rem_names_ids) + rem_names = np.char.add('InsertName', rem_names_ids) names = np.append(names, rem_names) for name in names: f.write(name + "\n") diff --git a/morphops/procrustes.py b/morphops/procrustes.py index 147379f..00f2cd4 100644 --- a/morphops/procrustes.py +++ b/morphops/procrustes.py @@ -20,7 +20,7 @@ import warnings def get_position(lmks): - """Returns the centroid of the set or sets of landmarks in `lmks`. + r"""Returns the centroid of the set or sets of landmarks in `lmks`. The centroid of a :math:`p` landmarks is simply the arithmetic mean of all the landmark positions. That is @@ -57,7 +57,7 @@ def get_position(lmks): return np.asarray(np.nanmean(lmks, axis=axis)) def get_scale(lmks): - """Returns the euclidean norm of the real matrix or matrices in `lmks`. + r"""Returns the euclidean norm of the real matrix or matrices in `lmks`. The euclidean norm of the real (p x k) matrix :math:`X` is calculated as @@ -182,7 +182,7 @@ def remove_scale(lmks, scale=None): return np.reshape(np.divide(lmks, scale_re), lmks_shape) def rotate(source, target, no_reflect=False): - """Rotates the landmark set `source` so as to minimize its sum of squared interlandmark distances to `target`. + r"""Rotates the landmark set `source` so as to minimize its sum of squared interlandmark distances to `target`. Say X=`source` and Y=`target`. By default :func:`rotate` tries to find @@ -260,25 +260,25 @@ def rotate(source, target, no_reflect=False): return result def opa(source, target, do_scaling=False, no_reflect=False): - """Performs Ordinary Procrustes Alignment to transform the landmark set + r"""Performs Ordinary Procrustes Alignment to transform the landmark set `source` such that the squared Euclidean distance between `source` and `target` is minimized. Say X=`source` and Y=`target` and `do_scaling` = `True`. :func:`opa` tries to find - .. math:: \operatorname*{argmin}_{\\beta > 0,\ R \in O(k),\ \gamma \in \mathbb{R}^k } D^2_{\mathtt{OPA}}(X, Y) = \| Y - \\beta X R - \mathbf{1_k} \gamma^T \|^2 + .. math:: \operatorname*{argmin}_{\beta > 0,\ R \in O(k),\ \gamma \in \mathbb{R}^k } D^2_{\mathtt{OPA}}(X, Y) = \| Y - \beta X R - \mathbf{1_k} \gamma^T \|^2 - If `do_scaling` = `False`, :math:`\\beta = 1`. If `no_reflect` = `True`, + If `do_scaling` = `False`, :math:`\beta = 1`. If `no_reflect` = `True`, then just as in :func:`rotate`, :func:`opa` will force :math:`R \in SO(k)`. The Ordinary (Procrustes) Sum of Squares or OSS is defined as - .. math:: OSS(X, Y) = \operatorname*{min}_{\\beta > 0,\ R \in O(k),\ \gamma \in \mathbb{R}^k } D^2_{\mathtt{OPA}}(X, Y) + .. math:: OSS(X, Y) = \operatorname*{min}_{\beta > 0,\ R \in O(k),\ \gamma \in \mathbb{R}^k } D^2_{\mathtt{OPA}}(X, Y) Note ---- - Generally for :func:`opa`, :math:`OSS(X1, X2) \\neq OSS(X2, X1)`. + Generally for :func:`opa`, :math:`OSS(X1, X2) \neq OSS(X2, X1)`. In contrast to :func:`opa`, :func:`gpa` is symmetric for the input matrices in that :math:`G(X1, X2) = G(X2, X1)`. @@ -297,8 +297,8 @@ def opa(source, target, do_scaling=False, no_reflect=False): do_scaling : bool, optional Flag indicating whether the best alignment should also find the optimal - :math:`\\beta` that minimizes :math:`D^2_{\mathtt{OPA}}`. The default - value of `do_scaling` is False, which means :math:`\\beta = 1`, or in + :math:`\beta` that minimizes :math:`D^2_{\mathtt{OPA}}`. The default + value of `do_scaling` is False, which means :math:`\beta = 1`, or in other words, `source` will not be scaled. no_reflect : bool, optional @@ -314,7 +314,7 @@ def opa(source, target, do_scaling=False, no_reflect=False): aligned to the `target`. b: numpy.float64 or int - A number representing the scaling factor :math:`\\beta` by which + A number representing the scaling factor :math:`\beta` by which `source` is scaled. R: numpy.ndarray @@ -328,7 +328,7 @@ def opa(source, target, do_scaling=False, no_reflect=False): oss: numpy.float64 This number represents the Ordinary (Procrustes) Sum of Squares, which is the minimum of :math:`D^2_{\mathtt{OPA}}`. Essentially, - the `oss` is the result of plugging in the optimal :math:`\\beta`, + the `oss` is the result of plugging in the optimal :math:`\beta`, :math:`R` and :math:`\gamma` into the :math:`D^2_{\mathtt{OPA}}` objective. @@ -396,17 +396,17 @@ def opa(source, target, do_scaling=False, no_reflect=False): def gpa(X, tol=1e-5,max_iters=10, do_project=False, do_scaling=False, no_reflect=False, unitize_mean=False): - """Performs Generalized Procrustes Alignment to transform all the landmark + r"""Performs Generalized Procrustes Alignment to transform all the landmark sets in `X` such that (a quantity proportional to) the sum of squared norms of pairwise differences between all the landmark sets is minimized. Say :code:`len(X) = n`. :func:`gpa` tries to find - .. math:: \operatorname*{argmin}_{\\beta_i > 0,\ R_i \in O(k),\ \gamma_i \in \mathbb{R}^k } g(X) = \\frac{1}{n} \sum_{i=1}^{n-1} { \sum_{j=i+1}^n {\| (\\beta_i X_i R_i + \mathbf{1_k} \gamma_i^T) - (\\beta_j X_j R_j + \mathbf{1_k} \gamma_j^T) \|^2}} + .. math:: \operatorname*{argmin}_{\beta_i > 0,\ R_i \in O(k),\ \gamma_i \in \mathbb{R}^k } g(X) = \frac{1}{n} \sum_{i=1}^{n-1} { \sum_{j=i+1}^n {\| (\beta_i X_i R_i + \mathbf{1_k} \gamma_i^T) - (\beta_j X_j R_j + \mathbf{1_k} \gamma_j^T) \|^2}} The Generalized (Procrustes) Sum of Squares or G is defined as - .. math:: G(X) = \operatorname*{inf}_{\\beta_i > 0,\ R_i \in O(k),\ \gamma_i \in \mathbb{R}^k } g(X) + .. math:: G(X) = \operatorname*{inf}_{\beta_i > 0,\ R_i \in O(k),\ \gamma_i \in \mathbb{R}^k } g(X) The GPA algorithm, per [drymar]_, tries to iteratively rotate and scale the landmark sets in `X` until the sum of squared differences is below `tol`. @@ -426,7 +426,7 @@ def gpa(X, tol=1e-5,max_iters=10, do_project=False, do_scaling=False, Note ---- - Generally for :func:`opa`, :math:`OSS(X1, X2) \\neq OSS(X2, X1)`. + Generally for :func:`opa`, :math:`OSS(X1, X2) \neq OSS(X2, X1)`. In contrast to :func:`opa`, :func:`gpa` is symmetric for the input matrices in that :math:`G(X1, X2) = G(X2, X1)`. @@ -453,8 +453,8 @@ def gpa(X, tol=1e-5,max_iters=10, do_project=False, do_scaling=False, has been achieved or `max_iters` is reached, whichever comes first. do_scaling : bool, optional - If `False`, :math:`\\beta_i = \\frac{1}{\| X'_i \|}`, where - :math:`X'_i` is the mean-centered :math:`X_i`. Else :math:`\\beta_i` is + If `False`, :math:`\beta_i = \frac{1}{\| X'_i \|}`, where + :math:`X'_i` is the mean-centered :math:`X_i`. Else :math:`\beta_i` is calculated as per [tenb]_. do_project: bool, optional @@ -482,12 +482,12 @@ def gpa(X, tol=1e-5,max_iters=10, do_project=False, do_scaling=False, b: numpy.ndarray A (n,)-shaped array representing the scaling factor - :math:`\\beta_i` by which the centered :math:`X'_i` is scaled. + :math:`\beta_i` by which the centered :math:`X'_i` is scaled. ssq: numpy.float64 This number represents the Generalized (Procrustes) Sum of Squares, which is the infinimum of :math:`g`. Essentially, - the `ssq` is the result of plugging in the optimal :math:`\\beta_i`, + the `ssq` is the result of plugging in the optimal :math:`\beta_i`, :math:`R_i` and :math:`\gamma_i` into the :math:`g` objective. Warns diff --git a/morphops/tests/helpers.py b/morphops/tests/helpers.py index d5e8cd2..c52d716 100644 --- a/morphops/tests/helpers.py +++ b/morphops/tests/helpers.py @@ -19,8 +19,8 @@ def make_haus(): The roof of the house is tilted to the left, which creates asymmetry. """ - haus = np.row_stack((make_ngon(4, 3*piOvr4), - make_ngon(3, -piOvr4/2) + [0, 1])) + haus = np.vstack((make_ngon(4, 3*piOvr4), + make_ngon(3, -piOvr4/2) + [0, 1])) haus_c = np.array([0,3/7]) haus0 = haus - haus_c haus0_b = np.sqrt(np.sum(np.square(haus0))) diff --git a/morphops/tests/test_tps.py b/morphops/tests/test_tps.py index 5bc1072..caa38a8 100644 --- a/morphops/tests/test_tps.py +++ b/morphops/tests/test_tps.py @@ -34,7 +34,7 @@ "zeros and A is a row of ones stacked on identity." ), (fourgon, fourgon_45, - np.zeros((4,2)),np.row_stack(([0,0],get_2d_rot(-np.pi/4))), + np.zeros((4,2)),np.vstack(([0,0],get_2d_rot(-np.pi/4))), "X is a square of side sqrt(2), Y = XR. Result should be that " "W is all zeros and A is a row of zeros stacked on R.T." )] diff --git a/morphops/tps.py b/morphops/tps.py index 693db96..8c23bb1 100644 --- a/morphops/tps.py +++ b/morphops/tps.py @@ -22,7 +22,7 @@ import warnings def K_matrix(X, Y=None): - """Calculates the upper-right (p,p) submatrix of the (p+k+1,p+k+1)-shaped + r"""Calculates the upper-right (p,p) submatrix of the (p+k+1,p+k+1)-shaped L matrix. Parameters @@ -45,11 +45,11 @@ def K_matrix(X, Y=None): In particular, if k = 2, then :math:`U(r) = r^2 \log(r^2)`, else :math:`U(r) = r`. - Note: Using :math:`\\alpha U(r)` instead of :math:`U(r)` for some - :math:`\\alpha \in \mathbb{R}` will not change the calculated spline. + Note: Using :math:`\alpha U(r)` instead of :math:`U(r)` for some + :math:`\alpha \in \mathbb{R}` will not change the calculated spline. Simple block matrix inverse formulae show that when calculating :math:`L^{-1}` - for the spline using :math:`\\alpha U(r)`, the non-uniform coefficients - multiplied to the :math:`U` terms will be scaled by :math:`\\frac{1}{\\alpha}` + for the spline using :math:`\alpha U(r)`, the non-uniform coefficients + multiplied to the :math:`U` terms will be scaled by :math:`\frac{1}{\alpha}` while the uniform coefficients will stay the same. """ num_coords = lmk_util.num_coords(X) @@ -154,7 +154,7 @@ def tps_coefs(X, Y): """ n_coords = lmk_util.num_coords(X) n_lmks = lmk_util.num_lmks(X) - Y_0 = np.row_stack((Y, np.zeros((n_coords+1,n_coords)))) + Y_0 = np.vstack((Y, np.zeros((n_coords+1,n_coords)))) L = L_matrix(X) Q = np.linalg.solve(L, Y_0) if np.any(np.isnan(Q)):