Non-exhaustive#

Non-exhaustive template aims to identify a global optimum without evaluating all possible configurations, e.g. via gradient-descent. However, whether non-exhaustive template matching is able to find a global optimum is highly dependent on the topology of the score space and available constraints. Typically, non-exhaustive template matching is useful when the rough orientation is known.

Within pytme, non-exhaustive template matching problems are distinguished based on whether they operate on numpy.ndarray objects or coordinates. Therefore, we categorize them into three different classes of matching problems:

Analogous to exhaustive template matching operations in pytme, the classes below identify the configuration maximizing the similarity between two densities.

FLC(target, template[, template_mask, ...])

Computes a normalized cross-correlation score of a target f a template g and a mask m:

The classes below maximize the similarity between a coordinate set with associated weights and a density. This case is commonly found when matching atomic structures to electron density maps.

CrossCorrelation(target, ...[, ...])

Computes the Cross-Correlation score as:

LaplaceCrossCorrelation(**kwargs)

Uses the same formalism as CrossCorrelation but with Laplace filtered weights (\(\nabla^{2}\)):

NormalizedCrossCorrelation(target, ...[, ...])

Computes a normalized version of the CrossCorrelation score based on the dot product of target_weights and template_weights, in order to reduce bias to regions of high local energy.

NormalizedCrossCorrelationMean(**kwargs)

Computes a similar score than NormalizedCrossCorrelation, but additionally factors in the mean of template and target.

MaskedCrossCorrelation(target, ...[, ...])

The Masked Cross-Correlation computes the similarity between target_weights and template_weights under respective masks.

PartialLeastSquareDifference(target, ...[, ...])

The Partial Least Square Difference (PLSQ) between the target \(f\) and the template \(g\) is calculated as:

Chamfer(target_coordinates, ...[, ...])

The Chamfer distance between the target \(f\) and the template \(g\) is calculated as:

MutualInformation(target, ...[, ...])

The Mutual Information (MI) score between the target \(f\) and the template \(g\) is calculated as:

Envelope([target_threshold])

The Envelope score (ENV) between the target \(f\) and the template \(g\) is calculated as:

NormalVectorScore(target_coordinates, ...[, ...])

The Normal Vector Score (NVS) between the target's \(f\) and the template \(g\)'s normal vectors is calculated as:

The classes below maximize the similarity between two coordinate sets and associated weights. This could be special points within densities, regular coordinate sets, or vectors between template and target that should be aligned.

Chamfer(target_coordinates, ...[, ...])

The Chamfer distance between the target \(f\) and the template \(g\) is calculated as:

NormalVectorScore(target_coordinates, ...[, ...])

The Normal Vector Score (NVS) between the target's \(f\) and the template \(g\)'s normal vectors is calculated as:

In order to perform non-exhaustive template matching, we have to initialize one of the objects outlined above. The following will showcase CrossCorrelation, but the procedure can be extended analogously to all scores shown above. We can use __call__ method of the created score_object to compute the score of the current configuration.

import numpy as np
from tme.matching_utils import create_mask
from tme.matching_optimization import CrossCorrelation

target = create_mask(
   mask_type="ellipse",
   radius=(15,10,5),
   shape=(51,51,51),
   center=(25,25,25),
)

# Note this could also be any other template that can be expressed as
# coordinates, e.g. an atomic structure.
template = create_mask(
   mask_type="ellipse",
   radius=(5,10,15),
   shape=(51,51,51),
   center=(32,17,25),
)
template_coordinates = np.array(np.where(template > 0))
template_weights = template[tuple(template_coordinates)]

score_object = CrossCorrelation(
   target=target,
   template_coordinates=template_coordinates,
   template_weights=template_weights,
   negate_score=True # Multiply returned score with -1 for minimization
)

score_object() # -543.0

Alternatively, score_object could have been instantiated using the wrapper function create_score_object().

create_score_object(score, **kwargs)

Initialize score object with name score using **kwargs.

pytme defines optimize_match(), which provides routines for non-exhaustive template matching given a score_object and bounds on possible translations and rotations.

optimize_match(score_object[, ...])

Find the translation and rotation optimizing the score returned by score_object with respect to provided bounds.

Given our previously defined score_object, we can find the optimal configuration:

from tme.matching_optimization import optimize_match

translation, rotation, refined_score = optimize_match(
   score_object=score_object,
   optimization_method="basinhopping",
   maxiter=50
)
translation
# array([-7.04120465,  8.17709819,  0.07964882])
rotation
# array([[ 3.0803247e-03,  1.7386347e-04,  9.9999523e-01],
#        [ 3.9309423e-02,  9.9922705e-01, -2.9481627e-04],
#        [-9.9922234e-01,  3.9310146e-02,  3.0711093e-03]],
#  dtype=float32)
refined_score #-3069.0

Note

optimize_match() will update score_object so using __call__ afterwards will return refined_score. Smaller values of maxiter might yield suboptimal results.

The computed translation and rotation can then be used to recover the configuration of the template with maximal similarity to the target. This is functionally similar to the operation performed using Density.match_densities and Density.match_structure_to_density.

from tme import Density

template_density = Density(template.astype(np.float32))
template_density = template_density.rigid_transform(
   translation = translation,
   rotation_matrix = rotation,
   order = 3,
   use_geometric_center=True
)

optimize_match() offers various optimization schemes to balance accuracy and runtime efficiency. Here we used optimization_method="basinhopping", which yields a good balance between those factors. For performance-critical applications, optimization_method="minimize" in conjunction with CrossCorrelation.grad can be used, which incorporates the analytical gradient of the scoring function into the optimization process. This approach can significantly speed up convergence, especially when the initial configuration is close to the global optimum.

The example below highlights this key characteristic. The initial orientation is too far from the correct one to be recovered without better initial values or tighter bounds on translation and rotation. Consequently, the optimization scheme becomes trapped in a local optimum. In practice, multiple initial conditions and potentially bounds should be evaluated to identify the global optimum.

score_object = CrossCorrelation(
   target=target,
   template_coordinates=template_coordinates,
   template_weights=template_weights,
   negate_score=True,
   return_gradient=True
)

translation, rotation, refined_score = optimize_match(
   score_object=score_object,
   optimization_method="minimize",
   maxiter=500,
   x0=[0,0,0,0,50,0]
)
translation
# array([-4.51057473,  9.10106869,  1.40713668])
rotation
# array([[ 0.09346937,  0.00730731,  0.99559534],
#       [-0.08880167,  0.9960488 ,  0.00102632],
#       [-0.99165404, -0.08850646,  0.09374895]],
#  dtype=float32)
refined_score #-2254.263671875

Tip

For custom optimization strategies, developers can utilize the CrossCorrelation.score() method, which accepts a tuple of translations and Euler angles as input.