"""
Similarity measures for the registration methods and factory to create similarity measures.
"""
from __future__ import print_function
from __future__ import absolute_import
from builtins import range
from builtins import object
from abc import ABCMeta, abstractmethod
import torch
from .data_wrapper import AdaptVal,USE_CUDA
from .data_wrapper import MyTensor
from . import utils
from math import floor
from .similarity_helper_omt import *
import torch.nn.functional as F
import numpy as np
from future.utils import with_metaclass
[docs]class SimilarityMeasure(with_metaclass(ABCMeta, object)):
"""Abstract base class for a similarity measure.
"""
def __init__(self, spacing, params):
self.spacing = spacing
"""pixel/voxel spacing"""
self.volumeElement = self.spacing.prod()
"""volume element"""
self.dim = len(spacing)
"""image dimension"""
self.params = params
"""external parameters"""
self.sigma = params['similarity_measure'][('sigma', 0.1, '1/sigma^2 is the weight in front of the similarity measure')]
"""1/sigma^2 is a balancing constant"""
[docs] def compute_similarity_multiNC(self, I0, I1, I0Source=None, phi=None):
"""
Compute the multi-image multi-channel image similarity between two images of format BxCxXxYzZ
:param I0: first image (the warped source image)
:param I1: second image (target image)
:param I0Source: source image (will typically not be used)
:param phi: map in the target image to warp the source image (will typically not be used)
:return: returns similarity measure
"""
sz = I0.size()
sim = torch.zeros(1).type_as(I0)
if I0Source is None and phi is None:
raise ValueError("no implemented in multi-batch-channel way1")
for nrI in range(sz[0]): # loop over all the images
sim = sim + self.compute_similarity_multiC(I0[nrI, ...], I1[nrI, ...], None, None )
elif I0Source is not None and phi is None:
# for nrI in range(sz[0]): # loop over all the images
# sim = sim + self.compute_similarity_multiC(I0[nrI, ...], I1[nrI, ...], I0Source[nrI,...], None)
sim = sim+self.compute_similarity(I0, I1, I0Source, None)
elif I0Source is None and phi is not None:
raise ValueError("no implemented in multi-batch-channel way3")
for nrI in range(sz[0]): # loop over all the images
sim = sim + self.compute_similarity_multiC(I0[nrI, ...], I1[nrI, ...], None, phi[nrI,...])
else:
# for nrI in range(sz[0]): # loop over all the images
# sim = sim + self.compute_similarity_multiC(I0[nrI, ...], I1[nrI, ...], I0Source[nrI,...], phi[nrI,...])
sim = sim+self.compute_similarity(I0, I1, I0Source, phi)
return sim
[docs] def compute_similarity_multiC(self, I0, I1, I0Source=None, phi=None):
"""
Compute the multi-channel image similarity between two images of format CxXxYzZ
:param I0: first image (the warped source image)
:param I1: second image (target image)
:param I0Source: source image (will typically not be used)
:param phi: map in the target image to warp the source image (will typically not be used)
:return: returns similarity measure
"""
sz = I0.size()
sim = torch.zeros(1).type_as(I0)
if I0Source is None:
for nrC in range(sz[0]): # loop over all the channels, just advect them all the same; if available map is the same for all channels
sim = sim + self.compute_similarity(I0[nrC, ...], I1[nrC, ...], None, phi)
else:
for nrC in range(sz[0]): # loop over all the channels, just advect them all the same; if available map is the same for all channels
sim = sim + self.compute_similarity(I0[nrC, ...], I1[nrC, ...],I0Source[nrC, ...],phi)
return sim/sz[0] # needs to be normalized based on the number of channels
[docs] @abstractmethod
def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Abstract method to compute the *single*-channel image similarity between two images of format XxYzZ.
This is the only method that should be overwritten by a specific implemented similarity measure.
The multi-channel variants then come for free.
For proper implementation it is important that the similarity measure is muliplied by
1./(self.sigma ** 2) and also by self.volumeElement if it is a volume integral
(and not a correlation measure for example)
:param I0: first image (the warped source image)
:param I1: second image (target image)
:param I0Source: source image (will typically not be used)
:param phi: map in the target image to warp the source image (will typically not be used)
:return: returns similarity measure
"""
pass
[docs] def set_sigma(self, sigma):
"""
Set balancing constant :math:`\\sigma`
:param sigma: balancing constant
"""
self.sigma = sigma
self.params['similarity_measure']['sigma']=sigma
[docs] def get_sigma(self):
"""
Get balancing constant
:return: balancing constant
"""
return self.sigma
[docs]class SSDSimilarity(SimilarityMeasure):
"""
Sum of squared differences (SSD) similarity measure.
:math:`1/sigma^2||I_0-I_1||^2`
"""
def __init__(self, spacing, params):
super(SSDSimilarity,self).__init__(spacing,params)
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the SSD measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
:return: SSD/sigma^2
"""
# TODO: This is to avoid a current pytorch bug 0.3.0 which cannot properly deal with infinity or NaN
return AdaptVal(utils.remove_infs_from_variable((I0 - I1) ** 2).sum() / (self.sigma ** 2) * self.volumeElement)
#return AdaptVal(((I0 - I1) ** 2).sum() / (self.sigma ** 2) * self.volumeElement)
[docs]class OptimalMassTransportSimilarity(SimilarityMeasure):
"""
Constructor
:param spacing: spacing of the grid (as in parameters structure)
:param std_dev: similarity measure parameter
:param sinkhorn_iterations: number of iterations in sinkhorn algorithm
:param std_sinkhorn: standard deviation of the entropic regularization (not too small to avoid nan)
"""
def __init__(self, spacing, params,sinkhorn_iterations = 200,std_sinkhorn = 0.08):
super(OptimalMassTransportSimilarity, self).__init__(spacing, params)
self.spacing = spacing
#self.params = params
self.std_dev = self.sigma
self.std_sinkhorn = std_sinkhorn
self.sinkhorn_iterations = sinkhorn_iterations
self.spline_order = params[('spline_order', 1, 'Spline interpolation order; 1 is linear interpolation (default); 3 is cubic spline')]
"""order of spline for interpolation (if needed)"""
[docs] def compute_similarity(self, I0, I1, I0Source, phi):
"""
Computes the SSD measure between two images
:param I0: first image (not used)
:param I1: second image (target image)
:param I0Source: source image (not warped)
:param phi: map to warp the source image to the target
:return: OMTSimilarity/sigma^2
"""
# todo: using the map here is not very efficient. It is much more efficient to apply the map
# todo: to an entire batch of images.
# FX: put your OMT code here; this is just a placeholder for now which is simple SSD (but using the source image and the map)
if phi is None:
raise ValueError('OptimalMassTransportSimiliary can only be computed for map-based models.')
# warp the source image (would be more efficient if we process a batch of images at once;
# but okay for now and no overhead if you only use one image pair at a time)
I1_warped = utils.compute_warped_image(I0Source, phi, self.spacing,self.spline_order)
# Encapsulate the data in tensor Variables
multiplier0 = torch.zeros(I0.size())
multiplier1 = torch.zeros(I1.size())
nr_iterations_sinkhorn = torch.Tensor([self.sinkhorn_iterations])
std_sink = torch.Tensor([self.std_sinkhorn])
# Compute the actual similarity
result = OTSimilarityHelper.apply(phi,I1_warped,I1,multiplier0,multiplier1,torch.Tensor(self.spacing),nr_iterations_sinkhorn,std_sink)
return result/(self.std_dev**2)
[docs]class NCCSimilarity(SimilarityMeasure):
"""
Computes a normalized-cross correlation based similarity measure between two images.
:math:`sim = (1-ncc^2)/(\\sigma^2)`
"""
def __init__(self, spacing, params):
super(NCCSimilarity,self).__init__(spacing,params)
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the NCC-based image similarity measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
:return: (1-NCC^2)/sigma^2
"""
# TODO: may require a safeguard against infinity
# this way of computing avoids the square root of the standard deviation
# I0mean = I0.mean()
# I1mean = I1.mean()
# I0_m_mean = I0-I0mean
# I1_m_mean = I1-I1mean
# nccSqr = (((I0_m_mean)*(I1_m_mean)).mean()**2)/\
# (((I0_m_mean)**2).mean()*((I1_m_mean)**2).mean())
n_batch = I0.shape[0]
dim = len(I0.shape[2:])
input_shape = [I0.shape[0], I0.shape[1], -1]+[1]*dim
I0 = I0.view(*input_shape)
I1 = I1.view(*input_shape)
I0mean = I0.mean(2)
I1mean = I1.mean(2)
I0_m_mean = I0-I0mean
I1_m_mean = I1-I1mean
nccSqr = (((I0_m_mean)*(I1_m_mean)).mean()**2)/\
(((I0_m_mean)**2).mean()*((I1_m_mean)**2).mean())
nccSqr =nccSqr.sum()
return AdaptVal((n_batch*1.-nccSqr)/self.sigma**2)
#ncc = ((I0-I0.mean().expand_as(I0))*(I1-I1.mean().expand_as(I1))).mean()/(I0.std()*I1.std())
# does not need to be multiplied by self.volumeElement (as we are dealing with a correlation measure)
#return AdaptVal((1.0-ncc**2) / (self.sigma ** 2))
[docs]class NCCPositiveSimilarity(SimilarityMeasure):
"""
Computes a normalized-cross correlation based similarity measure between two images. Only allows positive correlations.
:math:`sim = (1-ncc)/(\\sigma^2)`
"""
def __init__(self, spacing, params):
super(NCCPositiveSimilarity,self).__init__(spacing,params)
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the NCC-based image similarity measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
:return: (1-NCC)/sigma^2
"""
# TODO: may require a safeguard against infinity
# this way of computing avoids the square root of the standard deviation
n_batch = I0.shape[0]
dim = len(I0.shape[2:])
input_shape = [I0.shape[0], I0.shape[1], -1]+[1]*dim
I0 = I0.view(*input_shape)
I1 = I1.view(*input_shape)
I0mean = I0.mean(2)
I1mean = I1.mean(2)
I0_m_mean = I0 - I0mean
I1_m_mean = I1 - I1mean
ncc = (((I0_m_mean)*(I1_m_mean)).mean())/\
(torch.sqrt(((I0_m_mean)**2).mean())*torch.sqrt(((I1_m_mean)**2).mean()))
ncc = ncc.sum()
return AdaptVal((n_batch*1.-ncc)/self.sigma**2)
[docs]class NCCNegativeSimilarity(SimilarityMeasure):
"""
Computes a normalized-cross correlation based similarity measure between two images. Only allows negative correlations.
:math:`sim = (ncc)/(\\sigma^2)`
"""
def __init__(self, spacing, params):
super(NCCNegativeSimilarity,self).__init__(spacing,params)
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the NCC-based image similarity measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
:return: (NCC)/sigma^2
"""
# TODO: may require a safeguard against infinity
# this way of computing avoids the square root of the standard deviation
dim = len(I0.shape[2:])
input_shape = [I0.shape[0], I0.shape[1], -1] + [1] * dim
I0 = I0.view(*input_shape)
I1 = I1.view(*input_shape)
I0mean = I0.mean(2)
I1mean = I1.mean(2)
I0_m_mean = I0 - I0mean
I1_m_mean = I1 - I1mean
ncc = (((I0_m_mean) * (I1_m_mean)).mean()) / \
(torch.sqrt((I0_m_mean ** 2).mean()) * torch.sqrt((I1_m_mean ** 2).mean()))
ncc = ncc.sum()
return AdaptVal((ncc)/self.sigma**2)
[docs]class LNCCSimilarity(SimilarityMeasure):
"""This is an generalized LNCC; we implement multi-scale (means resolution)
multi kernel (means size of neighborhood) LNCC.
:param: resol_bound : type list, resol_bound[0]> resol_bound[1] >... resol_bound[end]
:param: kernel_size_ratio: type list, the ratio of the current input size
:param: kernel_weight_ratio: type list, the weight ratio of each kernel size, should sum to 1
:param: stride: type_list, the stride between each pixel that would compute its lncc
:param: dilation: type_list
Settings in json::
"similarity_measure": {
"develop_mod_on": false,
"sigma": 0.5,
"type": "lncc",
"lncc":{
"resol_bound":[-1],
"kernel_size_ratio":[[0.25]],
"kernel_weight_ratio":[[1.0]],
"stride":[0.25,0.25,0.25],
"dilation":[1]
}
For multi-scale multi kernel, e.g.,::
"resol_bound":[64,32],
"kernel_size_ratio":[[0.0625,0.125, 0.25], [0.25,0.5], [0.5]],
"kernel_weight_ratio":[[0.1,0.3,0.6],[0.3,0.7],[1.0]],
"stride":[0.25,0.25,0.25],
"dilation":[1,2,2] #[2,1,1]
or for single-scale single kernel, e.g.,::
"resol_bound":[-1],
"kernel_size_ratio":[[0.25]],
"kernel_weight_ratio":[[1.0]],
"stride":[0.25],
"dilation":[1]
Multi-scale is controlled by "resol_bound", e.g resol_bound = [128, 64], it means if input size>128, then it would compute multi-kernel
lncc designed for large image size, if 64<input_size<128, then it would compute multi-kernel lncc desiged for mid-size image, otherwise,
it would compute the multi-kernel lncc designed for small image.
Attention! we call it multi-scale just because it is designed for multi-scale registration or segmentation problem.
ONLY ONE scale would be activated during computing the similarity, which depends on the current input size.
At each scale, corresponding multi-kernel lncc is implemented, here multi-kernel means lncc with different window sizes
Loss = w1*lncc_win1 + w2*lncc_win2 ... + wn*lncc_winn, where /sum(wi) =1
for example. when (image size) S>128, three windows sizes can be used, namely S/16, S/8, S/4.
for easy notation, we use img_ratio to refer window size, the example here use the parameter [1./16,1./8,1.4]
In implementation, we compute lncc by calling convolution function, so in this case, the [S/16, S/8, S/4] refers
to the kernel size of convolution function. Intuitively, we would have another two parameters,
stride and dilation. For each window size (W), we recommend using W/4 as stride. In extreme case the stride can be 1, but
can large increase computation. The dilation expand the reception field, set dilation as 2 would physically twice the window size.
"""
def __init__(self, spacing, params):
super(LNCCSimilarity,self).__init__(spacing,params)
self.dim = len(spacing)
self.resol_bound = params['similarity_measure']['lncc'][('resol_bound',[128,64], "resolution bound for using different strategy")]
self.kernel_size_ratio = params['similarity_measure']['lncc'][('kernel_size_ratio',[[1./16,1./8,1./4],[1./4,1./2],[1./2]], "kernel size, ratio of input size")]
self.kernel_weight_ratio = params['similarity_measure']['lncc'][('kernel_weight_ratio',[[0.1, 0.3, 0.6],[0.3,0.7],[1.]], "kernel size, ratio of input size")]
self.strides = params['similarity_measure']['lncc'][('stride',[[1./4,1./4,1./4],[1./4,1./4],[1./4]], "step size, responded with ratio of kernel size")]
self.dilations = params['similarity_measure']['lncc'][('dilation',[[2,2,2],[2,2],[1]], "dilation param, responded with ratio of kernel size")]
if self.resol_bound[0] >-1:
assert len(self.resol_bound)+1 == len(self.kernel_size_ratio)
assert len(self.resol_bound)+1 == len(self.kernel_weight_ratio)
assert len(self.resol_bound)+1 == len(self.strides)
assert len(self.resol_bound)+1 == len(self.dilations)
def __stepup(self,img_sz):
max_scale = min(img_sz)
for i, bound in enumerate(self.resol_bound):
if max_scale >= bound:
self.kernel = [int(max_scale*kz) for kz in self.kernel_size_ratio[i]]
self.weight = self.kernel_weight_ratio[i]
self.stride = self.strides[i]
self.dilation = self.dilations[i]
break
if max_scale < self.resol_bound[-1]:
self.kernel = [int(max_scale*kz) for kz in self.kernel_size_ratio[-1]]
self.weight = self.kernel_weight_ratio[-1]
self.stride = self.strides[-1]
self.dilation = self.dilations[-1]
self.num_scale = len(self.kernel)
self.kernel_sz = [[k for _ in range(self.dim)] for k in self.kernel]
self.step = [[max(int((ksz + 1) * self.stride[scale_id]),1) for ksz in self.kernel_sz[scale_id]] for scale_id in range(self.num_scale)]
if USE_CUDA:
self.filter = [torch.ones([1, 1] + self.kernel_sz[scale_id]).cuda() for scale_id in range(self.num_scale)]
else:
self.filter = [torch.ones([1, 1] + self.kernel_sz[scale_id]) for scale_id in range(self.num_scale)]
if self.dim==1:
self.conv= F.conv1d
elif self.dim ==2:
self.conv= F.conv2d
elif self.dim ==3:
self.conv = F.conv3d
else:
raise ValueError(" Only 1-3d support")
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the NCC-based image similarity measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
"""
n_batch = I0.shape[0]
input = I0 #.view([1,1]+ list(I0.shape))
target =I1 #.view([1,1]+ list(I1.shape))
self.__stepup(img_sz=list(I0.shape[2:]))
input_2 = input ** 2
target_2 = target ** 2
input_target = input * target
lncc_total = 0.
for scale_id in range(self.num_scale):
input_local_sum = self.conv(input, self.filter[scale_id], padding=0, dilation=self.dilation[scale_id],
stride=self.step[scale_id]).view(input.shape[0], -1)
target_local_sum = self.conv(target, self.filter[scale_id], padding=0, dilation=self.dilation[scale_id],
stride=self.step[scale_id]).view(input.shape[0],
-1)
input_2_local_sum = self.conv(input_2, self.filter[scale_id], padding=0, dilation=self.dilation[scale_id],
stride=self.step[scale_id]).view(input.shape[0],
-1)
target_2_local_sum = self.conv(target_2, self.filter[scale_id], padding=0, dilation=self.dilation[scale_id],
stride=self.step[scale_id]).view(
input.shape[0], -1)
input_target_local_sum = self.conv(input_target, self.filter[scale_id], padding=0,
dilation=self.dilation[scale_id], stride=self.step[scale_id]).view(
input.shape[0], -1)
input_local_sum = input_local_sum.contiguous()
target_local_sum = target_local_sum.contiguous()
input_2_local_sum = input_2_local_sum.contiguous()
target_2_local_sum = target_2_local_sum.contiguous()
input_target_local_sum = input_target_local_sum.contiguous()
numel = float(np.array(self.kernel_sz[scale_id]).prod())
input_local_mean = input_local_sum / numel
target_local_mean = target_local_sum / numel
cross = input_target_local_sum - target_local_mean * input_local_sum - \
input_local_mean * target_local_sum + target_local_mean * input_local_mean * numel
input_local_var = input_2_local_sum - 2 * input_local_mean * input_local_sum + input_local_mean ** 2 * numel
target_local_var = target_2_local_sum - 2 * target_local_mean * target_local_sum + target_local_mean ** 2 * numel
lncc = cross * cross / (input_local_var * target_local_var + 1e-5)
lncc = 1 - lncc.mean()
lncc = n_batch* lncc
lncc_total += lncc * self.weight[scale_id]
return lncc_total / (self.sigma ** 2)
[docs]class LocalizedNCCSimilarity(SimilarityMeasure):
"""
Computes a normalized-cross correlation based similarity measure between two images.
:math:`sim = (1-ncc^2)/(\\sigma^2)`
"""
def __init__(self, spacing, params):
super(LocalizedNCCSimilarity,self).__init__(spacing,params)
#todo: maybe add some form of Gaussian weighing and tie it to the real image dimensions
self.gaussian_std = params['similarity_measure'][('gaussian_std', 0.025, 'standard deviation of Gaussian that will be used for local NCC computations')]
"""half the side length of the cube over which lNCC is computed"""
self.nr_of_elements_in_direction = None
self.weighting_coefficients = None
self.mask = None
self._create_gaussian_weighting(self.gaussian_std)
def _get_shifted_1d(self, I, x):
ret = torch.zeros_like(I)
sz = ret.size()
if x >= 0:
ret[0:sz[0] - x] = I[x:]
else: # x < 0
ret[-x:] = I[0:sz[0] + x]
return ret
def _get_shifted_2d(self,I,x,y):
ret = torch.zeros_like(I)
sz = ret.size()
if x>=0 and y>=0:
ret[0:sz[0]-x,0:sz[1]-y]=I[x:,y:]
elif x<0 and y>=0:
ret[-x:,0:sz[1]-y]=I[0:sz[0]+x,y:]
elif x>=0 and y<0:
ret[0:sz[0]-x,-y:] = I[x:, 0:sz[1]+y]
else: # x<0 and y<0
ret[-x:,-y:]= I[0:sz[0]+x,0:sz[1]+y]
return ret
def _get_shifted_3d(self, I, x, y, z):
ret = torch.zeros_like(I)
sz = ret.size()
if x >= 0 and y >= 0 and z>=0:
ret[0:sz[0] - x, 0:sz[1] - y,0:sz[2]-z] = I[x:, y:, z:]
elif x < 0 and y >= 0 and z>=0:
ret[-x:, 0:sz[1] - y,0:sz[2]-z] = I[0:sz[0] + x, y:, z:]
elif x >= 0 and y < 0 and z>=0:
ret[0:sz[0] - x, -y:,0:sz[2]-z] = I[x:, 0:sz[1] + y, z:]
elif x<0 and y<0 and z>=0:
ret[-x:, -y:,0:sz[2]-z] = I[0:sz[0] + x, 0:sz[1] + y, z:]
elif x >= 0 and y >= 0 and z<0:
ret[0:sz[0] - x, 0:sz[1] - y, -z:] = I[x:, y:, 0:sz[2]+z]
elif x < 0 and y >= 0 and z<0:
ret[-x:, 0:sz[1] - y, -z:] = I[0:sz[0] + x, y:, 0:sz[2]+z]
elif x >= 0 and y < 0 and z<0:
ret[0:sz[0] - x, -y:, -z:] = I[x:, 0:sz[1] + y, 0:sz[2]+z]
else: # x < 0 and y < 0 and z < 0:
ret[-x:, -y:, -z:] = I[0:sz[0] + x, 0:sz[1] + y, 0:sz[2]+z]
return ret
def _create_gaussian_weighting(self,sigma):
radiusSqr = (3.*sigma)**2
self.nr_of_elements_in_direction = MyTensor(self.dim).zero_()
for i in range(self.dim):
self.nr_of_elements_in_direction[i] = 3*sigma/self.spacing[i]
self.nr_of_elements_in_direction = torch.ceil(self.nr_of_elements_in_direction).int()
# now create the precomputed weights
self.weighting_coefficients = MyTensor(*list((2*self.nr_of_elements_in_direction+1).int())).zero_()
self.mask = torch.zeros_like(self.weighting_coefficients)
if self.dim==1:
for x in range(-self.nr_of_elements_in_direction[0],self.nr_of_elements_in_direction[0]+1):
currentRSqr = (x*self.spacing[0])**2
if currentRSqr<= radiusSqr:
self.weighting_coefficients[x+self.nr_of_elements_in_direction[0]]=np.exp(-currentRSqr/(2*sigma**2))
self.mask[x+self.nr_of_elements_in_direction[0]] = 1
# now normalize it
self.weighting_coefficients /= self.weighting_coefficients.sum()
elif self.dim==2:
for x in range(-self.nr_of_elements_in_direction[0],self.nr_of_elements_in_direction[0]+1):
for y in range(-self.nr_of_elements_in_direction[1],self.nr_of_elements_in_direction[1]+1):
currentRSqr = (x*self.spacing[0])**2 + (y*self.spacing[1])**2
if currentRSqr<= radiusSqr:
self.weighting_coefficients[x+self.nr_of_elements_in_direction[0],y+self.nr_of_elements_in_direction[1]]=\
np.exp(-currentRSqr/(2*sigma**2))
self.mask[x+self.nr_of_elements_in_direction[0],y+self.nr_of_elements_in_direction[1]] = 1
# now normalize it
self.weighting_coefficients /= self.weighting_coefficients.sum()
elif self.dim==3:
for x in range(-self.nr_of_elements_in_direction[0],self.nr_of_elements_in_direction[0]+1):
for y in range(-self.nr_of_elements_in_direction[1],self.nr_of_elements_in_direction[1]+1):
for z in range(-self.nr_of_elements_in_direction[2], self.nr_of_elements_in_direction[2] + 1):
currentRSqr = (x*self.spacing[0])**2 + (y*self.spacing[1])**2 + (z*self.spacing[2])**2
if currentRSqr<= radiusSqr:
self.weighting_coefficients[x+self.nr_of_elements_in_direction[0],
y+self.nr_of_elements_in_direction[1],
z+self.nr_of_elements_in_direction[2]]= np.exp(-currentRSqr/(2*sigma**2))
self.mask[x+self.nr_of_elements_in_direction[0],
y+self.nr_of_elements_in_direction[1],
z+self.nr_of_elements_in_direction[2]] = 1
# now normalize it
self.weighting_coefficients /= self.weighting_coefficients.sum()
else:
raise ValueError('Only dimensions 1,2, and 3 are supported.')
def _compute_local_squared_cross_correlation(self,I0,I1):
ones = torch.ones_like(I0)
sumOnes = torch.zeros_like(I0)
sumI0 = torch.zeros_like(I0)
sumI1 = torch.zeros_like(I0)
sumI0I1 = torch.zeros_like(I0)
sumI0I0 = torch.zeros_like(I0)
sumI1I1 = torch.zeros_like(I0)
I0I0 = I0*I0
I1I1 = I1*I1
I0I1 = I0*I1
if self.dim==1:
for x in range(-self.nr_of_elements_in_direction[0],self.nr_of_elements_in_direction[0]+1):
if self.mask[self.nr_of_elements_in_direction[0] + x] > 0:
current_weight = self.weighting_coefficients[self.nr_of_elements_in_direction[0]+x]
sumOnes += current_weight*self._get_shifted_1d(ones, x)
sumI0 += current_weight*self._get_shifted_1d(I0, x)
sumI1 += current_weight*self._get_shifted_1d(I1, x)
sumI0I1 += current_weight*self._get_shifted_1d(I0I1, x)
sumI0I0 += current_weight*self._get_shifted_1d(I0I0, x)
sumI1I1 += current_weight*self._get_shifted_1d(I1I1, x)
elif self.dim==2:
for x in range(-self.nr_of_elements_in_direction[0],self.nr_of_elements_in_direction[0]+1):
for y in range(-self.nr_of_elements_in_direction[1],self.nr_of_elements_in_direction[1]+1):
if self.mask[self.nr_of_elements_in_direction[0]+x,self.nr_of_elements_in_direction[1]+y]>0:
current_weight = self.weighting_coefficients[self.nr_of_elements_in_direction[0]+x,self.nr_of_elements_in_direction[1]+y]
sumOnes += current_weight*self._get_shifted_2d(ones,x,y)
sumI0 += current_weight*self._get_shifted_2d(I0,x,y)
sumI1 += current_weight*self._get_shifted_2d(I1,x,y)
sumI0I1 += current_weight*self._get_shifted_2d(I0I1,x,y)
sumI0I0 += current_weight*self._get_shifted_2d(I0I0,x,y)
sumI1I1 += current_weight*self._get_shifted_2d(I1I1,x,y)
elif self.dim==3:
for x in range(-self.nr_of_elements_in_direction[0], self.nr_of_elements_in_direction[0] + 1):
for y in range(-self.nr_of_elements_in_direction[1], self.nr_of_elements_in_direction[1] + 1):
for z in range(-self.nr_of_elements_in_direction[2], self.nr_of_elements_in_direction[2] + 1):
if self.mask[
self.nr_of_elements_in_direction[0] + x,
self.nr_of_elements_in_direction[1] + y,
self.nr_of_elements_in_direction[2] + z] > 0:
current_weight = self.weighting_coefficients[
self.nr_of_elements_in_direction[0] + x,
self.nr_of_elements_in_direction[1] + y,
self.nr_of_elements_in_direction[2] + z]
sumOnes += current_weight*self._get_shifted_3d(ones, x, y, z)
sumI0 += current_weight*self._get_shifted_3d(I0, x, y, z)
sumI1 += current_weight*self._get_shifted_3d(I1, x, y, z)
sumI0I1 += current_weight*self._get_shifted_3d(I0I1, x, y, z)
sumI0I0 += current_weight*self._get_shifted_3d(I0I0, x, y, z)
sumI1I1 += current_weight*self._get_shifted_3d(I1I1, x, y, z)
else:
raise ValueError('Only supported in dimensions 1, 2, and 3')
# 1/n\sum_i (I0-mean(I0))(I1-mean(I1)) = 1/n \sum_i (I0I1 -I0 mean(I1) - mean(I0)I1 + mean(I0)mean(I1) )
# ... = ( 1/n \sum_i I0 I1 ) - mean(I0)mean(I1)
# \sigma_0 = 1/n \sum_i (I0-mean(I0))^2 = 1/n \sum_i I0^2 - 2I0 mean(I0) + mean(I0)^2
# ... = (1/n \sum_i I0^2 ) - mean(I0)^2
meanI0 = sumI0/sumOnes
meanI1 = sumI1/sumOnes
nom = sumI0I1/sumOnes - meanI0*meanI1
sig0Sqr = (sumI0I0/sumOnes - meanI0**2)
sig1Sqr = (sumI1I1/sumOnes - meanI1**2)
# todo: maybe find a little less hacky solution to deal with division by zero
# we are returning the square here, because it is squared later anyway
# and taking the sub-gradient for the square root at zero is not well defined
eps = 1e-6 # to avoid division by zero
lnccSqr = (nom*nom+eps)/(sig0Sqr*sig1Sqr+eps)
return lnccSqr
[docs] def compute_similarity(self, I0, I1, I0Source=None, phi=None):
"""
Computes the NCC-based image similarity measure between two images
:param I0: first image
:param I1: second image
:param I0Source: not used
:param phi: not used
:return: (1-NCC^2)/sigma^2
"""
# TODO: may require a safeguard against infinity
lnccSqr = self._compute_local_squared_cross_correlation(I0,I1)
# does not need to be multiplied by self.volumeElement (as we are dealing with a correlation measure)
sim_measure = AdaptVal((1.0-lnccSqr).sum() / (I0.numel()*self.sigma ** 2))
#print( 'sim_measure = ' + str( sim_measure.data.numpy()))
return sim_measure
[docs]class SimilarityMeasureFactory(object):
"""
Factory to quickly generate similarity measures that can then be used by the different registration algorithms.
"""
def __init__(self,spacing):
self.spacing = spacing
"""image spacing"""
self.dim = len( spacing )
"""dimension of image"""
self.similarity_measure_default_type = 'ssd'
"""default image similarity measure"""
self.simMeasures = {
'ssd': SSDSimilarity,
'ncc': NCCSimilarity,
'ncc_positive': NCCPositiveSimilarity,
'ncc_negative': NCCNegativeSimilarity,
'lncc': LNCCSimilarity,#LocalizedNCCSimilarity,
'omt': OptimalMassTransportSimilarity
}
"""currently implemented similiarity measures"""
[docs] def add_similarity_measure(self,simName,simClass):
"""
Adds a new custom similarity measure
:param simName: desired name of the similarity measure
:param simClass: similiarity measure class (whcih can be instantiated)
"""
print('Registering new similarity measure ' + simName + ' in factory')
self.simMeasures[simName]=simClass
[docs] def print_available_similarity_measures(self):
"""
Prints all the available similarity measures
"""
print(self.simMeasures)
[docs] def set_similarity_measure_default_type_to_ssd(self):
"""
Set the default similarity measure to SSD
"""
self.similarity_measure_default_type = 'ssd'
[docs] def set_similarity_measure_default_type_to_omt(self):
"""
Set the default similarity measure to OMT (optimal mass transport)
"""
self.similarity_measure_default_type = 'omt'
[docs] def set_similarity_measure_default_type_to_ncc(self):
"""
Set the default similarity measure to NCC
"""
self.similarity_measure_default_type = 'ncc'
[docs] def set_similarity_measure_default_type_to_ncc_positive(self):
"""
Set the default similarity measure to positive NCC (i.e., only positive correlations allowed)
"""
self.similarity_measure_default_type = 'ncc_positive'
[docs] def set_similarity_measure_default_type_to_ncc_negative(self):
"""
Set the default similarity measure to positive NCC (i.e., only negative correlations allowed)
"""
self.similarity_measure_default_type = 'ncc_negative'
[docs] def set_similarity_measure_default_type_to_lncc(self):
"""
Set the default similarity measure to *localized* NCC
"""
self.similarity_measure_default_type = 'lncc'
[docs] def create_similarity_measure(self, params):
"""
Create the actual similarity measure
:param params: ParameterDict() object holding the parameters which can contol similarity measure settings
:return: returns a similarity measure (which can then be used to evaluate similarities)
"""
cparams = params[('similarity_measure',{},'settings for the similarity measure')]
similarityMeasureType = cparams[('type', self.similarity_measure_default_type, 'type of similarity measure (ssd/ncc)')]
if similarityMeasureType in self.simMeasures:
print('Using ' + similarityMeasureType + ' similarity measure')
return self.simMeasures[similarityMeasureType](self.spacing,params)
else:
raise ValueError( 'Similarity measure: ' + similarityMeasureType + ' not known')