参考：
https://blog.mbedded.ninja/programming/signal-processing/image-processing/image-registration/modality-independent-neighbourhood-descriptor-mind/
《MIND: Modality independent neighbourhood descriptor for multi-modal deformable registration》论文
2D版的MIND_shchojj的博客-CSDN博客
MIND（模态独立邻域描述符）是一种图像配准算法，旨在提取每个像素/体素周围局部区域的结构信息。
一、原理介绍
1、公式介绍
MIND表达式如下：
Dp(I,x,x+r)指以x点为中心的Patch块和以x+r为中心的Patch块的距离，即两个Patch块之差的平方。
V(I,x)是指图像中x点的6邻域搜索空间内各点的Dp(I,x,x+r)相加再平均。(这个邻域搜索空间后面介绍)
MIND(I,x,r)是感兴趣点x的邻域搜索空间内这n个点的Dp(I,x,x+r)除以该点的V(I,x)，再负对数。
邻域搜索空间neigh Search Space:分为三种
紧密型采样（全采样） 稀疏性采样（45°采样） 6邻域采样
*红色小块为感兴趣点，即需要求MIND特征的像素点。
如果设置一个5*5的一个neigh,紧密型采样的邻域搜索空间包括除红色像素以外的所有灰色像素点；稀疏型采样同理如中图所示；6邻域采样应该是文章针对三维图像来讲的，对于二维的图像，应该也可以说是4邻域。
Patch:3*3的patch是指一个3*3的像素小块
2、图示说明
从一个5*5的像素块开始理解：
*蓝色像素4的地方是感兴趣点，即正要计算MIND特征的像素点，
*靛青色像素是设置的4邻域搜索空间
开始计算之前，先设置参数patch(3*3),neigh(3*3,4邻域采样)
（1）4邻域搜索空间内各点的Dp(I,x,x+r)
向量x=(1,2)，向量x+r={(1,1),(0,2),(1,3),(2,2)}
同理求完四个点的Dp(I,x,x+r)如下图所示：
（2）求V(I,x)
V(I,x)在2维图像中需要求4邻域搜索空间的4个点的Dp(I,x,x+r)的平均
到这里只求得了一个像素点（1,2）的V(I,x)，对于整个5*5的像素块或者256*256的图像，我们就需要求出这个像素块或者图像的V(I,x)的一个矩阵。
array([[ 15.33, 41.97, 68.19, 87.61, 60.97],[ 47.06, 93.78, 126.25, 137.28, 90.56],[ 93.72, 141.75, 157.06, 129.78, 81.75],[109.42, 147.36, 152.08, 99.03, 61.08],[ 77.69, 95.56, 94.03, 49.36, 31.5 ]])（3）计算MIND特征
那么MIND(I,(1,2))=[0.300,0.275,0.361,0.616]
同样到这里，我们只获得了点（1，2）的MIND特征，而对于这个5*5的像素块，我们还要获得其余24个点的MIND特征。
最终得到的这个5*5的像素块的MIND特征如下：
mind_descriptors =array([[[0.199, 0.404, 0.258, 0.884],[0.142, 0.554, 0.381, 0.611],[0.354, 0.301, 0.408, 0.421],[0.505, 0.446, 0.265, 0.307],[0.955, 0.374, 0.205, 0.25 ]],[[0.367, 0.615, 0.126, 0.643],[0.176, 0.604, 0.266, 0.649],[0.3 , 0.275, 0.361, 0.616],[0.395, 0.33 , 0.328, 0.428],[0.878, 0.244, 0.248, 0.344]],[[0.43 , 0.845, 0.142, 0.354],[0.301, 0.572, 0.255, 0.416],[0.326, 0.339, 0.377, 0.44 ],[0.419, 0.257, 0.552, 0.308],[0.765, 0.251, 0.445, 0.214]],[[0.49 , 0.887, 0.224, 0.188],[0.376, 0.589, 0.308, 0.269],[0.349, 0.388, 0.371, 0.365],[0.383, 0.199, 0.525, 0.46 ],[0.621, 0.211, 0.413, 0.339]],[[0.488, 0.948, 0.325, 0.121],[0.518, 0.558, 0.39 , 0.162],[0.432, 0.512, 0.412, 0.201],[0.575, 0.202, 0.575, 0.274],[0.528, 0.42 , 0.459, 0.18 ]]])如果想要比较一个5*5和像素块A和一个5*5的像素块B的结构相似性，那么就要先获得A和B的形状为（5,5,4）的MIND特征矩阵。然后将每个像素点的对应的MIND矩阵加和作为该像素点的MIND特征，最后就又会得到A的5*5的MIND特征图和B的MIND特征图。论文中的MIND特征图展示如下：
二、Python实现代码Pytorch（2D版）
import torch
import numpy as np
import torchvision
import SimpleITK as sitk
import torch.nn.functional as F
import matplotlib.pyplot as plt
def gaussian_kernel(sigma, sz):xpos_vec = np.arange(sz)ypos_vec = np.arange(sz)output = np.ones([1, 1,sz, sz], dtype=np.single)midpos = sz // 2for xpos in xpos_vec:for ypos in ypos_vec:output[:,:,xpos,ypos] = np.exp(-((xpos-midpos)**2 + (ypos-midpos)**2) / (2 * sigma**2)) / (2 * np.pi * sigma**2)return output
def torch_image_translate(input_, tx, ty, interpolation='nearest'):# got these parameters from solving the equations for pixel translations# on https://www.tensorflow.org/api_docs/python/tf/contrib/image/transformtranslation_matrix = torch.zeros([1, 3, 3], dtype=torch.float)translation_matrix[:, 0, 0] = 1.0translation_matrix[:, 1, 1] = 1.0translation_matrix[:, 0, 2] = -2*tx/(input_.size()[2]-1)translation_matrix[:, 1, 2] = -2*ty/(input_.size()[3]-1)translation_matrix[:, 2, 2] = 1.0grid = F.affine_grid(translation_matrix[:, 0:2, :], input_.size()).to(input_.device)wrp = F.grid_sample(input_, grid, mode=interpolation)return wrp
def Dp(image, xshift, yshift, sigma, patch_size):shift_image = torch_image_translate(image, xshift, yshift, interpolation='nearest')#将image在x、y方向移动I`(a)diff = torch.sub(image, shift_image)#计算差分图I-I`(a)diff_square = torch.mul(diff, diff)#(I-I`(a))^2res = torch.conv2d(diff_square, weight =torch.from_numpy(gaussian_kernel(sigma, patch_size)), stride=1, padding=3)#C*(I-I`(a))^2return resdef MIND(image, patch_size, neigh_size, sigma, eps,image_size0,image_size1, name='MIND'):# compute the Modality independent neighbourhood descriptor (MIND) of input image.# suppose the neighbor size is R, patch size is P.# input image is 384 x 256 x input_c_dim# output MIND is (384-P-R+2) x (256-P-R+2) x R*Rreduce_size = int((patch_size + neigh_size - 2) / 2)#卷积后减少的size# estimate the local variance of each pixel within the input image.Vimg = torch.add(Dp(image, -1, 0, sigma, patch_size), Dp(image, 1, 0, sigma, patch_size))Vimg = torch.add(Vimg, Dp(image, 0, -1, sigma, patch_size))Vimg = torch.add(Vimg, Dp(image, 0, 1, sigma, patch_size))#sum(Dp)Vimg = torch.div(Vimg,4) + torch.mul(torch.ones_like(Vimg), eps)#防除零# estimate the (R*R)-length MIND feature by shifting the input image by R*R times.xshift_vec = np.arange( -(neigh_size//2), neigh_size - (neigh_size//2))#邻域计算yshift_vec = np.arange(-(neigh_size // 2), neigh_size - (neigh_size // 2))#邻域计算iter_pos = 0for xshift in xshift_vec:for yshift in yshift_vec:if (xshift,yshift) == (0,0):continueMIND_tmp = torch.exp(torch.mul(torch.div(Dp(image, xshift, yshift, sigma, patch_size), Vimg), -1))#exp(-D(I)/V(I))tmp = MIND_tmp[:, :, reduce_size:(image_size0 - reduce_size), reduce_size:(image_size1 - reduce_size)]if iter_pos == 0:output = tmpelse:output = torch.cat([output,tmp], 1)iter_pos = iter_pos + 1# normalization.input_max, input_indexes = torch.max(output, dim=1)output = torch.div(output,input_max)return output
def abs_criterion(in_, target):return torch.mean(torch.abs(in_ - target))
if __name__ == '__main__':patch_size=7neigh_size=9sigma=2.0eps=1e-5image_size0=512image_size1=512ct_path = r'F:\dataset\coarseReg\00001_SRO1_ct_axial_drr.nii.gz'mr_path = r'F:\dataset\coarseReg\00001_SRO1_mr_axial_drr.nii.gz'ct_sitk = sitk.ReadImage(ct_path,sitk.sitkFloat32)mr_sitk = sitk.ReadImage(mr_path,sitk.sitkFloat32)ct_axial_drr = sitk.GetArrayFromImage(ct_sitk)[np.newaxis ,np.newaxis,:,:]mr_axial_drr = sitk.GetArrayFromImage(mr_sitk)[np.newaxis ,np.newaxis,:,: ]plt.imshow(ct_axial_drr.squeeze())plt.show()ct_axial_drr = torch.from_numpy(ct_axial_drr)mr_axial_drr = torch.from_numpy(mr_axial_drr)A_MIND = MIND(ct_axial_drr, patch_size, neigh_size, sigma, eps,image_size0,image_size1, name='realA_MIND')B_MIND = MIND(mr_axial_drr, patch_size, neigh_size, sigma, eps,image_size0,image_size1, name='realA_MIND')g_loss_MIND = abs_criterion(A_MIND, B_MIND)print('g_loss_MIND', g_loss_MIND)# tf = np.load( r"C:\Users\shcho\Desktop\validation\tf.npy" )# tc = np.load( r"C:\Users\shcho\Desktop\validation\tc.npy" )# plt.imshow(tf)# plt.show()# plt.imshow(tc)# plt.show()# print(tf.shape, tc.shape)## # print(tf.transpose((2,0,1)).shape, tc.shape)# t = tf.transpose((2,0,1)) -tc## # t = tf -tc## print(np.max(t))可以通过上面代码提取整张图像的MIND特征，且该函数可以被用于基于CycleGAN网络处理MRtoCT图像的损失函数并得到很好的结果，如论文《Unsupervised MR-to-CT Synthesis using Structure Constrained CycleGAN》所述，同理可用于相似的CycleGAN网络模型处理医学图像的任务中来约束A，B域的结构信息。
