python中的加权相关函数

Question

我试图根据链接中的文章，公式2，实现加权相关函数：

http://staff.ustc.edu.cn/~lshao/papers/paper07.pdf

假设n个元素中有3个向量s、r和w。 矢量w由以下公式得出：

w = |r|/(1+D)
D = |s - k*r|
k = (r_Transpose * s)/(r_Transpose*r)

我想实现文章中描述的加权相关函数的公式。我的执行是否正确？ 我从维度[224640]的矩阵开始，这意味着我有224个元素的640个向量。我想计算这640个向量与另一个向量r之间的加权相关系数。这640个向量中的每一个都是向量s

ref = reference  
ref_mean = np.mean(ref)  # calcolo il valore medio dello spettro di riferimento
sens = 190  


frame_correlation = np.zeros((1,640))
img_correlation = np.zeros((nf,npixels))

for i in range(nf):
    frame_test = dati_new[:,:,i]   Selection of one matrix from a cell of matrices
    for j in range(npixels):

        spettro_test = frame_test[:,j]     # is my vector s
        spettro_test = np.reshape(spettro_test,(224,1))
        spettro_test_mean = np.mean(spettro_test) 

        k = np.dot(np.transpose(ref),spettro_test)/np.dot(np.transpose(ref),ref)
        k = k[0][0]
        D = np.abs(spettro_test - k*ref)
        W = np.abs(ref)/(1+D)

        # NUMERATOR OF FORMULA IN THE ARTICLE
        numeratore = np.sum(W*(spettro_test - spettro_test_mean)*(ref - ref_mean))

        # DENOMINATOR
        den1_ex = np.sqrt(np.sum(W*np.power(spettro_test - spettro_test_mean,2)))
        den2_ex = np.sqrt(np.sum(W*np.power(ref  - ref_mean,2)))
        denominatore = den1_ex * den2_ex
        rho = numeratore/denominatore

        if rho < 0:
            rho  = 0
        if rho  > 1: # for safety reason
            rho = 1

        if rho >=0.99:
            rho = (sens*rho)/100
        frame_correlation[:,j]= rho
    img_correlation[i,:] = frame_correlation

Answer 1

这是我编写的代码，用于实现从矩阵中选择的两个数组之间的加权相关函数

ref = reference  
ref_mean = np.mean(ref)  
sens = 190

nf = n #number of matrices    
frame_correlation = np.zeros((1,640))
img_correlation = np.zeros((nf,npixels))

for i in range(nf):
    frame_test = dati_new[:,:,i]  #dati_new is a 3D structure made of nf matrices    
    for j in range(npixels):

        spettro_test = frame_test[:,j]    
        spettro_test = np.reshape(spettro_test,(224,1))
        spettro_test_mean = np.mean(spettro_test) 

        # calcolo del peso per lo spettro selezionato
        k = np.dot(np.transpose(ref),spettro_test)/np.dot(np.transpose(ref),ref)
        k = k[0][0]
        D = np.abs(spettro_test - k*ref)
        W = np.abs(ref)/(1+D)

        # Definizione del numeratore del coefficiente di correlazione

        numeratore = np.sum(W*(spettro_test - spettro_test_mean)*(ref - ref_mean))

        # Definizione del denominatore del coefficiente di correlazione

        den1_ex = np.sqrt(np.sum(W*np.power(spettro_test - spettro_test_mean,2)))
        den2_ex = np.sqrt(np.sum(W*np.power(ref  - ref_mean,2)))
        denominatore = den1_ex * den2_ex
        rho = numeratore/denominatore

        if rho < 0:
            rho  = 0
        if rho  > 1: # just in case
            rho = 1

        if rho >=0.998:
            rho = (sens*rho)/100
        frame_correlation[:,j]= rho
    img_correlation[i,:] = frame_correlation

img_correlation = np.array(img_correlation)

fig, ax=plt.subplots()
ax.imshow(img_correlation,cmap="gray", origin="lower")
plt.title('correlation coefficient image')
plt.xlabel("Pixels")
plt.ylabel("Number of frames")
plt.show()