CODEWITHSUNDEEP
python
tang
tang

Reputation: 31

Calculate the mutual information value of 3D image

I would like to ask, how to use Python to calculate the mutual information value and standardized mutual information value of two three-dimensional images? Any help would be appreciated. I have tried to use the function in sklearn to calculate the mutual information value a, but I found that when the mutual information is calculated again after rotating or translating the image, the value B is obtained. The difference between a and B is very small, so I personally think this method is not applicable to images. Am I right?Here's the code I used with the sklearn package.

def NMI(img1,img2):
    img1_ = sitk.GetArrayFromImage(img1)
    img2_ = sitk.GetArrayFromImage(img2)
    img2_ = np.reshape(img2_, -1)
    img1_ = np.reshape(img1_, -1)
    normalized_mutual_infor = mr.normalized_mutual_info_score(img1_, img2_)
    nmi = normalized_mutual_infor 
    print(nmi)

fixed_image = sitk.ReadImage(r"D:\Lung CT\RIDER Lung CT\001_1.mha", sitk.sitkFloat32)
moving_image = sitk.ReadImage(r"D:\Lung CT\RIDER Lung CT\001_2.mha", sitk.sitkFloat32)
tfm1 = sitk.ReadTransform(r'D:\6freedom\1_text2.tfm')
x = tfm1.GetParameters()[3]
y = tfm1.GetParameters()[4]
z = tfm1.GetParameters()[5]

transform1 = sitk.Euler3DTransform(tfm1)
transform1.SetParameters((0, 0, 0, x, y, z))
resample = sitk.Resample(moving_image, fixed_image, transform1, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
NMI(fixed_image, resample)
#When the parameter is (0,0,0, x, y, z), the result is 0.524628297588729
#When the parameter is (1,0,0, x, y, z), the result is 0.4657578384754303
#The unit of rotation is radians, so the image has been rotated a lot, 
#but the difference between the two results is very small.

Upvotes: 0

Views: 1540

Answers (1)

Ibraheem
Ibraheem

Reputation: 55

From this nice notebook, it seems one can use the joint histogram of the input images e.g.

import numpy as np
def mutual_information(hgram):
     # Mutual information for joint histogram
     # Convert bins counts to probability values
     pxy = hgram / float(np.sum(hgram))
     px = np.sum(pxy, axis=1) # marginal for x over y
     py = np.sum(pxy, axis=0) # marginal for y over x
     px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
     # Now we can do the calculation using the pxy, px_py 2D arrays
     nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
     return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
hist_2d, x_edges, y_edges = np.histogram2d(img1.ravel(),img2.ravel(),bins=20)
mi = mutual_information(hist_2d)
print(mi)

Upvotes: 2

Related Questions