Page MenuHomec4science
Files F65887190

point_matching.py
No OneTemporary
Actions

File Metadata

Created
Thu, Jun 6, 21:23

point_matching.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
from scipy.optimize import linear_sum_assignment as LSA
from rrompy.utilities.base.types import Tuple, Np1D, Np2D, HFEng
__all__ = ['pointMatching', 'cutOffDistance', 'angleTable',
'chordalMetricTable', 'chordalMetricAdjusted']
def pointMatching(distanceMatrix:Np2D) -> Np1D:
return LSA(distanceMatrix)[1]
def cutOffDistance(x:Np1D, murange : Tuple[float, float] = [- 1., 1.]) -> Np1D:
mu0 = np.mean(murange)
musig = murange[0] - mu0
isInf = np.isinf(x)
dist = np.empty(len(x))
dist[isInf] = np.inf
xEffR = x[np.logical_not(isInf)] - mu0
if np.isclose(musig, 0.):
dist[np.logical_not(isInf)] = np.abs(xEffR)
else:
xEffR /= musig
bernEff = (xEffR ** 2. - 1) ** .5
dist[np.logical_not(isInf)] = np.max(np.vstack((
np.abs(xEffR + bernEff), np.abs(xEffR - bernEff)
)), axis = 0) - 1.
return dist
def angleTable(X:Np2D, Y:Np2D, HFEngine : HFEng = None,
is_state : bool = True) -> Np2D:
if HFEngine is None:
innerT = Y.dot(X.T.conj())
normX = np.linalg.norm(X, axis = 1)
normY = np.linalg.norm(Y, axis = 1)
else:
innerT = HFEngine.innerProduct(X, Y, is_state = is_state)
normX = HFEngine.norm(X, is_state = is_state)
normY = HFEngine.norm(Y, is_state = is_state)
xInf = np.where(np.isclose(normX, 0.))[0]
normX[xInf] = 1.
return - (np.abs(innerT / normX).T - normY)
def chordalMetricTable(x:Np1D, y:Np1D) -> Np2D:
x, y = np.array(x), np.array(y)
xInf, yInf = np.where(np.isinf(x))[0], np.where(np.isinf(y))[0]
x[xInf], y[yInf] = 0., 0.
T = np.abs(np.tile(x.reshape(-1, 1), len(y)) - y.reshape(1, -1))
T[xInf, :], T[:, yInf] = 1., 1.
T[np.ix_(xInf, yInf)] = 0.
T = (T.T * (np.abs(x) ** 2. + 1.) ** -.5).T * (np.abs(y) ** 2. + 1.) ** -.5
return T
def chordalMetricAdjusted(x:Np1D, y:Np1D, w : float = 0, X : Np2D = None,
Y : Np2D = None, HFEngine : HFEng = None,
is_state : bool = True) -> Np2D:
dist = chordalMetricTable(x, y)
if w == 0: return dist
distAdj = angleTable(X, Y, HFEngine, is_state)
return dist + w * distAdj

Event Timeline

c4science · Help