metrics

Loss functions for registration and reconstruction tasks

Image similarity metrics

Compute the similarity between a fixed X-ray \(\mathbf I\) and a moving X-ray \(\mathbf{\hat I}\), where \(\mathbf{\hat I}\) is rendered from an estimated camera pose (registration) or volume (reconstruction).

We implement patchwise variants of the following metrics:

  • Normalized Cross Correlation (NCC)
  • Multiscale Normalized Cross Correlation (mNCC)
  • Gradient Normalized Cross Correlation (gNCC)
Tip

If patch_size=None, the similarity metric is computed over the entire image.

source

NormalizedCrossCorrelation2d

def NormalizedCrossCorrelation2d(
    patch_size:NoneType=None, eps:float=1e-05
):

Compute Normalized Cross Correlation between two batches of images.

source

MultiscaleNormalizedCrossCorrelation2d

def MultiscaleNormalizedCrossCorrelation2d(
    patch_sizes:list=[None], patch_weights:list=[1.0], eps:float=1e-05
):

Compute Normalized Cross Correlation between two batches of images at multiple scales.

source

GradientNormalizedCrossCorrelation2d

def GradientNormalizedCrossCorrelation2d(
    patch_size:NoneType=None, sigma:float=1.0, **kwargs
):

Compute Normalized Cross Correlation between the image gradients of two batches of images.

source

MutualInformation

def MutualInformation(
    sigma:float=0.1, num_bins:int=256, epsilon:float=1e-10, normalize:bool=True
):

Mutual Information.

Geodesic distances for SE(3)

One can define geodesic pseudo-distances on \(\mathbf{SO}(3)\) and \(\mathbf{SE}(3)\).1 This let’s us measure registration error (in radians and millimeters, respectively) on poses, rather than needed to compute the projection of fiducials.

We implement two geodesics on \(\mathbf{SE}(3)\):

  • The logarithmic geodesic
  • The double geodesic

Logarithmic Geodesic

Given two rotation matrices \(\mathbf R_A, \mathbf R_B \in \mathbf{SO}(3)\), the angular distance between their axes of rotation is

\[ d_\theta(\mathbf R_A, \mathbf R_B) = \arccos \left( \frac{\mathrm{trace}(\mathbf R_A^T \mathbf R_B) - 1}{2} \right) = \| \log (\mathbf R_A^T \mathbf R_B) \| \,, \]

where \(\log(\cdot)\) is the logarithm map on \(\mathbf{SO}(3)\).2 Using the logarithm map on \(\mathbf{SE}(3)\), this generalizes to a geodesic loss function on camera poses \({\mathbf T}_A, {\mathbf T}_B \in \mathbf{SE}(3)\):

\[ \mathcal L_{\mathrm{log}}({\mathbf T}_A, {\mathbf T}_B) = \| \log({\mathbf T}_A^{-1} {\mathbf T}_B) \| \,. \]

source

LogGeodesicSE3

def LogGeodesicSE3():

Calculate the distance between transforms in the log-space of SE(3).

# SE(3) distance
geodesic_se3 = LogGeodesicSE3()

pose_1 = convert(
    torch.tensor([[0.1, 1.0, torch.pi]]),
    torch.ones(1, 3),
    parameterization="euler_angles",
    convention="ZYX",
)
pose_2 = convert(
    torch.tensor([[0.1, 1.1, torch.pi]]),
    torch.zeros(1, 3),
    parameterization="euler_angles",
    convention="ZYX",
)

geodesic_se3(pose_1, pose_2)
tensor([1.7354])

Double Geodesic

We can also formulate a geodesic distance on \(\mathbf{SE}(3)\) with units of length. Using the camera’s focal length \(f\), we convert the angular distance to an arc length:

\[ d_\theta(\mathbf R_A, \mathbf R_B; f) = \frac{f}{2} d_\theta(\mathbf R_A, \mathbf R_B) \,. \]

When combined with the Euclidean distance on the translations \(d_t(\mathbf t_A, \mathbf t_B) = \| \mathbf t_A - \mathbf t_B \|\), this yields the double geodesic loss on \(\mathbf{SE}(3)\):3

\[ \mathcal L_{\mathrm{geo}}({\mathbf T}_A, {\mathbf T}_B; f) = \sqrt{d^2_\theta(\mathbf R_A, \mathbf R_B; f) + d^2_t(\mathbf t_A, \mathbf t_B)} \,. \]

source

DoubleGeodesicSE3

def DoubleGeodesicSE3(
    sdd:float, # Source-to-detector distance
    eps:float=1e-06, # Avoid overflows in sqrt
):

Calculate the angular and translational geodesics between two SE(3) transformation matrices.

# Angular distance and translational distance both in mm
double_geodesic = DoubleGeodesicSE3(1020.0)

pose_1 = convert(
    torch.tensor([[0.1, 1.0, torch.pi]]),
    torch.ones(1, 3),
    parameterization="euler_angles",
    convention="ZYX",
)
pose_2 = convert(
    torch.tensor([[0.1, 1.1, torch.pi]]),
    torch.zeros(1, 3),
    parameterization="euler_angles",
    convention="ZYX",
)

double_geodesic(pose_1, pose_2)  # Angular, translational, double geodesics
(tensor([51.0000]), tensor([1.7321]), tensor([51.0294]))