gattia / cycpd Goto Github PK

from cycpd import affine_registration

def affine_cpd(X, Y, max_iterations=100, init_only=False):
    reg = affine_registration(**{'X': X, 'Y': Y, 'max_iterations': max_iterations, 'print_reg_params': False, 'verbose': False})
    # B is the rotation matrix, t is the translation vector
    if not init_only:
        B, t = reg.register()

    return reg

File "cycpd\cython\cython_functions.pyx", line 209, in cython_functions.__pyx_fused_cpdef
TypeError: No matching signature found

def _memoryview_safe(x):
    """Make array safe to run in a Cython memoryview-based kernel. These
    kernels typically break down with the error ``ValueError: buffer source
    array is read-only`` when running in dask distributed.
    """
    if not x.flags.writeable:
        if not x.flags.owndata:
            x = x.copy(order='C')
        x.setflags(write=True)
    return x

def register(self, callback=lambda **kwargs: None):
    tic = time.time()
    self.transform_point_cloud()
    self.X = _memoryview_safe(self.X)
    self.Y = _memoryview_safe(self.Y)
    self.sigma2 = cy.initialize_sigma2(self.X, self.Y)
    ......

def expectation(self):
    self.X = _memoryview_safe(self.X)
    self.TY = _memoryview_safe(self.TY)
    self.P1, self.Pt1, self.PX, self.Np, self.E = cy.expectation_2(self.X,
                                                                   self.TY,
                                                                   self.sigma2,
                                                                   self.M,
                                                                   self.N,
                                                                   self.D,
                                                                   self.w)

@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
@cython.cdivision(True) # Do division using C?
def expectation_3(my_type[:,:] X, my_type[:,:] TY, my_type sigma2, int M, int N, int D, double w):
    # TY is MxD (same as Y - it is just transformed Y). M = rows of points Y, and D is dimensions.
    # tic = time.time()
    # cdef Py_ssize_t x_i_shape, x_j_shape, y_i_shape, y_j_shape
    if my_type is double:
        dtype = np.double

    cdef Py_ssize_t n, m, d  # i, j, k

    cdef double ksig, diff, w_tmp, den, tmp_total, Np, E # den = sp in original matlab, tmp_total = razn
    # cdef double [:] temp_x
    temp_x = np.zeros(D, dtype=np.double)
    cdef double [:] temp_x_view = temp_x

    ksig = -2.0 * sigma2
    w_tmp = (w * M * c_pow(-ksig * 3.14159265358979,0.5*D))/((1-w)*N) # c in the original manuscript

    cdef my_type[:, :] X_view = X
    cdef my_type[:, :] TY_view = TY
	
    probP = np.zeros((M, N), dtype=dtype)
    cdef my_type[:,:] probP_view = probP

     # Make P array
    # cdef my_type[:] P
    P = np.zeros(M, dtype=dtype)
    cdef my_type[:] P_view = P

    # Make P1 array
    # cdef my_type[:] P1
    P1= np.zeros(M, dtype=dtype)
    cdef my_type[:] P1_view = P1

    # Make Pt1 array
    # cdef my_type[:] Pt1
    Pt1= np.zeros(N, dtype=dtype)
    cdef my_type[:] Pt1_view = Pt1

    # Make Pt1 array
    # cdef my_type[:, :] Px
    Px= np.zeros((M, D), dtype=dtype)
    cdef my_type[:, :] Px_view = Px


    for n in range(N):
        den = 0
        for m in range(M):
            tmp_total = 0
            for d in range(D):
                diff = X_view[n, d] - TY_view[m, d]
                diff = diff * diff
                tmp_total += diff
            P_view[m] = c_exp(tmp_total/ksig)
            probP_view[m, n] = c_exp(tmp_total/ksig)
            den += P_view[m]
        den += w_tmp
        for m in range(M):
            probP_view[m, n] = probP_view[m, n] / den
        Pt1_view[n] = 1-w_tmp/den

        for d in range(D):
            temp_x_view[d] = X_view[n, d] / den

        for m in range(M):
            P1_view[m] += P_view[m] / den

            for d in range(D):
                Px_view[m, d] += temp_x_view[d] * P_view[m]
        E += -c_log(den)

    Np = 0
    for m in range(M):
        Np += P1_view[m]

    E += D * Np * c_log(sigma2)/2

    return probP, P1, Pt1, Px, Np, E