时间戳

首页 > 解决方案 > python中的加权相关函数

python - python中的加权相关函数

问题描述

我尝试根据链接中的文章，公式编号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)

我想实现文章中描述的加权相关函数的公式。我的实施是否正确？我从维度为 [224,640] 的矩阵开始，这意味着我有 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

标签： pythonfor-loopcorrelation

解决方案

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

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()

推荐阅读