This documentation page refers to a previous release of DIALS (1.14).
Click here to go to the corresponding page for the latest version of DIALS

Source code for dials.algorithms.indexing.fft3d

#!/usr/bin/env python
# -*- mode: python; coding: utf-8; indent-tabs-mode: nil; python-indent: 2 -*-
#
# dials.algorithms.indexing.fft3d.py
#
#  Copyright (C) 2014 Diamond Light Source
#
#  Author: Richard Gildea
#
#  This code is distributed under the BSD license, a copy of which is
#  included in the root directory of this package.

from __future__ import absolute_import, division
from __future__ import print_function
import math

import libtbx
from dials.util import Sorry
from scitbx import fftpack
from scitbx import matrix
from cctbx import crystal, uctbx, xray
from cctbx.array_family import flex
from dials.algorithms.indexing.indexer import indexer_base, optimise_basis_vectors
from dials.algorithms.indexing.indexer import (
    vector_group,
    is_approximate_integer_multiple,
)
from dxtbx.model.experiment_list import Experiment, ExperimentList

import logging

logger = logging.getLogger(__name__)


[docs]class indexer_fft3d(indexer_base): def __init__(self, reflections, imagesets, params): super(indexer_fft3d, self).__init__(reflections, imagesets, params)
[docs] def find_lattices(self): if self.params.multiple_lattice_search.cluster_analysis_search: self.find_basis_vector_combinations_cluster_analysis() self.debug_show_candidate_basis_vectors() if self.params.debug_plots: self.debug_plot_candidate_basis_vectors() crystal_models = self.candidate_crystal_models if self.params.multiple_lattice_search.max_lattices is not None: crystal_models = crystal_models[ : self.params.multiple_lattice_search.max_lattices ] else: self.find_candidate_basis_vectors() self.debug_show_candidate_basis_vectors() if self.params.debug_plots: self.debug_plot_candidate_basis_vectors() self.candidate_crystal_models = self.find_candidate_orientation_matrices( self.candidate_basis_vectors ) crystal_model, n_indexed = self.choose_best_orientation_matrix( self.candidate_crystal_models ) if crystal_model is not None: crystal_models = [crystal_model] else: crystal_models = [] experiments = ExperimentList() for cm in crystal_models: for imageset in self.imagesets: experiments.append( Experiment( imageset=imageset, beam=imageset.get_beam(), detector=imageset.get_detector(), goniometer=imageset.get_goniometer(), scan=imageset.get_scan(), crystal=cm, ) ) return experiments
[docs] def map_centroids_to_reciprocal_space_grid(self): d_min = self.params.fft3d.reciprocal_space_grid.d_min n_points = self.gridding[0] rlgrid = 2 / (d_min * n_points) # real space FFT grid dimensions cell_lengths = [n_points * d_min / 2 for i in range(3)] self.fft_cell = uctbx.unit_cell(cell_lengths + [90] * 3) self.crystal_symmetry = crystal.symmetry( unit_cell=self.fft_cell, space_group_symbol="P1" ) logger.info("FFT gridding: (%i,%i,%i)" % self.gridding) grid = flex.double(flex.grid(self.gridding), 0) selection = self.reflections["id"] == -1 if self.params.b_iso is libtbx.Auto: self.params.b_iso = -4 * d_min ** 2 * math.log(0.05) logger.debug("Setting b_iso = %.1f" % self.params.b_iso) from dials.algorithms.indexing import map_centroids_to_reciprocal_space_grid map_centroids_to_reciprocal_space_grid( grid, self.reflections["rlp"], selection, d_min, b_iso=self.params.b_iso ) reflections_used_for_indexing = selection.iselection() self.reciprocal_space_grid = grid self.reflections_used_for_indexing = reflections_used_for_indexing
[docs] def fft(self): if self.params.fft3d.reciprocal_space_grid.d_min is libtbx.Auto: # rough calculation of suitable d_min based on max cell # see also Campbell, J. (1998). J. Appl. Cryst., 31(3), 407-413. # fft_cell should be greater than twice max_cell, so say: # fft_cell = 2.5 * max_cell # then: # fft_cell = n_points * d_min/2 # 2.5 * max_cell = n_points * d_min/2 # a little bit of rearrangement: # d_min = 5 * max_cell/n_points max_cell = self.params.max_cell d_min = 5 * max_cell / self.params.fft3d.reciprocal_space_grid.n_points d_spacings = 1 / self.reflections["rlp"].norms() self.params.fft3d.reciprocal_space_grid.d_min = max(d_min, min(d_spacings)) logger.info( "Setting d_min: %.2f" % self.params.fft3d.reciprocal_space_grid.d_min ) n_points = self.params.fft3d.reciprocal_space_grid.n_points self.gridding = fftpack.adjust_gridding_triple( (n_points, n_points, n_points), max_prime=5 ) n_points = self.gridding[0] self.map_centroids_to_reciprocal_space_grid() self.d_min = self.params.fft3d.reciprocal_space_grid.d_min logger.info( "Number of centroids used: %i" % ((self.reciprocal_space_grid > 0).count(True)) ) # gb_to_bytes = 1073741824 # bytes_to_gb = 1/gb_to_bytes # (128**3)*8*2*bytes_to_gb # 0.03125 # (256**3)*8*2*bytes_to_gb # 0.25 # (512**3)*8*2*bytes_to_gb # 2.0 fft = fftpack.complex_to_complex_3d(self.gridding) grid_complex = flex.complex_double( reals=self.reciprocal_space_grid, imags=flex.double(self.reciprocal_space_grid.size(), 0), ) grid_transformed = fft.forward(grid_complex) # self.grid_real = flex.pow2(flex.abs(grid_transformed)) self.grid_real = flex.pow2(flex.real(grid_transformed)) # self.grid_real = flex.pow2(flex.imag(self.grid_transformed)) del grid_transformed if self.params.debug: self.debug_write_ccp4_map(map_data=self.grid_real, file_name="fft3d.map") if self.params.fft3d.peak_search == "flood_fill": self.find_peaks() elif self.params.fft3d.peak_search == "clean": self.find_peaks_clean()
[docs] def find_peaks(self): grid_real_binary = self.grid_real.deep_copy() rmsd = math.sqrt( flex.mean( flex.pow2( grid_real_binary.as_1d() - flex.mean(grid_real_binary.as_1d()) ) ) ) grid_real_binary.set_selected( grid_real_binary < (self.params.rmsd_cutoff) * rmsd, 0 ) grid_real_binary.as_1d().set_selected(grid_real_binary.as_1d() > 0, 1) grid_real_binary = grid_real_binary.iround() from cctbx import masks flood_fill = masks.flood_fill(grid_real_binary, self.fft_cell) if flood_fill.n_voids() < 4: # Require at least peak at origin and one peak for each basis vector raise Sorry( "Indexing failed: fft3d peak search failed to find sufficient number of peaks." ) # the peak at the origin might have a significantly larger volume than the # rest so exclude any anomalously large peaks from determining minimum volume from scitbx.math import five_number_summary outliers = flex.bool(flood_fill.n_voids(), False) grid_points_per_void = flood_fill.grid_points_per_void() min_x, q1_x, med_x, q3_x, max_x = five_number_summary(grid_points_per_void) iqr_multiplier = 5 iqr_x = q3_x - q1_x cut_x = iqr_multiplier * iqr_x outliers.set_selected(grid_points_per_void.as_double() > (q3_x + cut_x), True) # print q3_x + cut_x, outliers.count(True) isel = ( grid_points_per_void > int( self.params.fft3d.peak_volume_cutoff * flex.max(grid_points_per_void.select(~outliers)) ) ).iselection() if self.params.optimise_initial_basis_vectors: self.volumes = flood_fill.grid_points_per_void().select(isel) sites_cart = flood_fill.centres_of_mass_cart().select(isel) sites_cart_optimised = optimise_basis_vectors( self.reflections["rlp"].select(self.reflections_used_for_indexing), sites_cart, ) self.sites = self.fft_cell.fractionalize(sites_cart_optimised) diffs = sites_cart_optimised - sites_cart norms = diffs.norms() flex.min_max_mean_double(norms).show() perm = flex.sort_permutation(norms, reverse=True) for p in perm[:10]: logger.debug(sites_cart[p], sites_cart_optimised[p], norms[p]) # only use those vectors which haven't shifted too far from starting point sel = norms < (5 * self.fft_cell.parameters()[0] / self.gridding[0]) self.sites = self.sites.select(sel) self.volumes = self.volumes.select(sel) # diff = (self.sites - flood_fill.centres_of_mass_frac().select(isel)) # flex.min_max_mean_double(diff.norms()).show() else: self.sites = flood_fill.centres_of_mass_frac().select(isel) self.volumes = flood_fill.grid_points_per_void().select(isel)
[docs] def find_candidate_basis_vectors(self): self.fft() # hijack the xray.structure class to facilitate calculation of distances xs = xray.structure(crystal_symmetry=self.crystal_symmetry) for i, site in enumerate(self.sites): xs.add_scatterer(xray.scatterer("C%i" % i, site=site)) xs = xs.sites_mod_short() sites_cart = xs.sites_cart() lengths = flex.double([matrix.col(sc).length() for sc in sites_cart]) perm = flex.sort_permutation(lengths) xs = xs.select(perm) volumes = self.volumes.select(perm) if self.params.debug: with open("peaks.pdb", "wb") as f: f.write(xs.as_pdb_file()) sites_frac = xs.sites_frac() vectors = xs.sites_cart() norms = vectors.norms() sel = (norms > self.params.min_cell) & (norms < (2 * self.params.max_cell)) vectors = vectors.select(sel) vectors = [matrix.col(v) for v in vectors] volumes = volumes.select(sel) # XXX loop over these vectors and sort into groups similar to further down # group similar angle and lengths, also catch integer multiples of vectors vector_groups = [] relative_length_tolerance = 0.1 angle_tolerance = 5 # degrees orth = self.fft_cell.orthogonalize for v, volume in zip(vectors, volumes): length = v.length() if length < self.params.min_cell or length > (2 * self.params.max_cell): continue matched_group = False for group in vector_groups: mean_v = group.mean() mean_v_length = mean_v.length() if ( abs(mean_v_length - length) / max(mean_v_length, length) < relative_length_tolerance ): angle = mean_v.angle(v, deg=True) if angle < angle_tolerance: group.append(v, length, volume) matched_group = True break elif abs(180 - angle) < angle_tolerance: group.append(-v, length, volume) matched_group = True break if not matched_group: group = vector_group() group.append(v, length, volume) vector_groups.append(group) vectors = [g.mean() for g in vector_groups] volumes = flex.double(max(g.volumes) for g in vector_groups) ## sort by length # lengths = flex.double([v.length() for v in vectors]) # perm = flex.sort_permutation(lengths) # sort by peak size perm = flex.sort_permutation(volumes, reverse=True) volumes = volumes.select(perm) vectors = [vectors[i] for i in perm] for i, (v, volume) in enumerate(zip(vectors, volumes)): logger.debug("%s %s %s" % (i, v.length(), volume)) lengths = flex.double(v.length() for v in vectors) perm = flex.sort_permutation(lengths) # exclude vectors that are (approximately) integer multiples of a shorter # vector unique_vectors = [] unique_volumes = flex.double() for p in perm: v = vectors[p] is_unique = True for i, v_u in enumerate(unique_vectors): if (unique_volumes[i] > volumes[p]) and is_approximate_integer_multiple( v_u, v ): logger.debug( "rejecting %s: integer multiple of %s" % (v.length(), v_u.length()) ) is_unique = False break if is_unique: unique_vectors.append(v) unique_volumes.append(volumes[p]) # re-sort by peak volume perm = flex.sort_permutation(unique_volumes, reverse=True) vectors = [unique_vectors[i] for i in perm] volumes = unique_volumes.select(perm) # for i, (v, volume) in enumerate(zip(vectors, volumes)): # logger.debug("%s %s %s" %(i, v.length(), volume)) self.candidate_basis_vectors = vectors return self.candidate_basis_vectors
[docs] def find_basis_vector_combinations_cluster_analysis(self): self.fft() # hijack the xray.structure class to facilitate calculation of distances xs = xray.structure(crystal_symmetry=self.crystal_symmetry) for i, site in enumerate(self.sites): xs.add_scatterer(xray.scatterer("C%i" % i, site=site)) xs = xs.sites_mod_short() xs = xs.select(xs.sites_frac().norms() < 0.45) cell_multiplier = 10 xs1 = xs.customized_copy( unit_cell=uctbx.unit_cell( [xs.unit_cell().parameters()[0] * cell_multiplier] * 3 ) ) xs1.set_sites_cart(xs.sites_cart()) xs = xs1 sites_cart = xs.sites_cart() lengths = flex.double([matrix.col(sc).length() for sc in sites_cart]) xs = xs.select(flex.sort_permutation(lengths)) if self.params.debug: with open("peaks.pdb", "wb") as f: f.write(xs.as_pdb_file()) vector_heights = flex.double() sites_frac = xs.sites_frac() pair_asu_table = xs.pair_asu_table(distance_cutoff=self.params.max_cell) asu_mappings = pair_asu_table.asu_mappings() distances = crystal.calculate_distances(pair_asu_table, sites_frac) vectors = [] difference_vectors = [] pairs = [] for di in distances: if di.distance < self.params.min_cell: continue i_seq, j_seq = di.i_seq, di.j_seq if i_seq > j_seq: continue pairs.append((i_seq, j_seq)) rt_mx_ji = di.rt_mx_ji site_frac_ji = rt_mx_ji * sites_frac[j_seq] site_cart_ji = xs.unit_cell().orthogonalize(site_frac_ji) site_cart_i = xs.unit_cell().orthogonalize(sites_frac[i_seq]) vectors.append(matrix.col(site_cart_ji)) diff_vec = matrix.col(site_cart_i) - matrix.col(site_cart_ji) if diff_vec[0] < 0: # only one hemisphere of difference vector space diff_vec = -diff_vec difference_vectors.append(diff_vec) params = self.params.multiple_lattice_search.cluster_analysis if params.method == "dbscan": i_cluster = self.cluster_analysis_dbscan(difference_vectors) min_cluster_size = 1 elif params.method == "hcluster": i_cluster = self.cluster_analysis_hcluster(difference_vectors) i_cluster -= 1 # hcluster starts counting at 1 min_cluster_size = params.min_cluster_size if self.params.debug_plots: self.debug_plot_clusters( difference_vectors, i_cluster, min_cluster_size=min_cluster_size ) clusters = [] min_cluster_size = params.min_cluster_size for i in range(max(i_cluster) + 1): isel = (i_cluster == i).iselection() if len(isel) < min_cluster_size: continue clusters.append(isel) cluster_point_sets = [] centroids = [] cluster_sizes = flex.int() difference_vectors = flex.vec3_double(difference_vectors) from libtbx.utils import flat_list for cluster in clusters: points = flat_list([pairs[i] for i in cluster]) cluster_point_sets.append(set(points)) d_vectors = difference_vectors.select(cluster) cluster_sizes.append(len(d_vectors)) centroids.append(d_vectors.mean()) # build a graph where each node is a centroid from the difference vector # cluster analysis above, and an edge is defined when there is a # significant overlap between the sets of peaks in the FFT map that # contributed to the difference vectors in two clusters import networkx as nx G = nx.Graph() G.add_nodes_from(range(len(cluster_point_sets))) cutoff_frac = 0.25 for i in range(len(cluster_point_sets)): for j in range(i + 1, len(cluster_point_sets)): intersection_ij = cluster_point_sets[i].intersection( cluster_point_sets[j] ) union_ij = cluster_point_sets[i].union(cluster_point_sets[j]) frac_connected = len(intersection_ij) / len(union_ij) if frac_connected > cutoff_frac: G.add_edge(i, j) # iteratively find the maximum cliques in the graph # break from the loop if there are no cliques remaining or there are # fewer than 3 vectors in the remaining maximum clique # Allow 1 basis vector to be shared between two cliques, to allow for # cases where two lattices share one basis vectors (e.g. two plate # crystals exactly aligned in one direction, but not in the other two) distinct_cliques = [] cliques = list(nx.find_cliques(G)) cliques = sorted(cliques, key=len, reverse=True) for i, clique in enumerate(cliques): clique = set(clique) if len(clique) < 3: break is_distinct = True for c in distinct_cliques: if len(c.intersection(clique)) > 1: is_distinct = False break if is_distinct: distinct_cliques.append(clique) this_set = set() for i_cluster in clique: this_set = this_set.union(cluster_point_sets[i_cluster]) logger.info("Clique %i: %i lattice points" % (i + 1, len(this_set))) assert len(distinct_cliques) > 0 logger.info("Estimated number of lattices: %i" % len(distinct_cliques)) self.candidate_basis_vectors = [] self.candidate_crystal_models = [] for clique in distinct_cliques: sel = flex.size_t(list(clique)) vectors = flex.vec3_double(centroids).select(sel) perm = flex.sort_permutation(vectors.norms()) vectors = [matrix.col(vectors[p]) for p in perm] # exclude vectors that are (approximately) integer multiples of a shorter # vector unique_vectors = [] for v in vectors: is_unique = True for v_u in unique_vectors: if is_approximate_integer_multiple( v_u, v, relative_tolerance=0.01, angular_tolerance=0.5 ): is_unique = False break if is_unique: unique_vectors.append(v) vectors = unique_vectors self.candidate_basis_vectors.extend(vectors) candidate_orientation_matrices = self.find_candidate_orientation_matrices( vectors ) if len(candidate_orientation_matrices) == 0: continue crystal_model, n_indexed = self.choose_best_orientation_matrix( candidate_orientation_matrices ) if crystal_model is None: continue # map to minimum reduced cell crystal_symmetry = crystal.symmetry( unit_cell=crystal_model.get_unit_cell(), space_group=crystal_model.get_space_group(), ) cb_op = crystal_symmetry.change_of_basis_op_to_minimum_cell() crystal_model = crystal_model.change_basis(cb_op) self.candidate_crystal_models.append(crystal_model) if self.params.debug: file_name = "vectors.pdb" a = self.params.max_cell cs = crystal.symmetry(unit_cell=(a, a, a, 90, 90, 90), space_group="P1") xs = xray.structure(crystal_symmetry=cs) for v in difference_vectors: v = matrix.col(v) xs.add_scatterer(xray.scatterer("C", site=v / (a / 10))) xs.sites_mod_short() with open(file_name, "wb") as f: f.write(xs.as_pdb_file()) for crystal_model in self.candidate_crystal_models: logger.debug(crystal_model)
[docs] def cluster_analysis_hcluster(self, vectors): from hcluster import linkage, fcluster import numpy params = self.params.multiple_lattice_search.cluster_analysis.hcluster X = numpy.array(vectors) linkage_method = params.linkage.method linkage_metric = params.linkage.metric criterion = params.cutoff_criterion Z = linkage(X, method=linkage_method, metric=linkage_metric) cutoff = params.cutoff i_cluster = fcluster(Z, cutoff, criterion=criterion) i_cluster = flex.int(i_cluster.astype(numpy.int32)) return i_cluster
[docs] def cluster_analysis_dbscan(self, vectors): import numpy as np from sklearn.cluster import DBSCAN from sklearn.preprocessing import StandardScaler vectors = flex.vec3_double(vectors) X = np.array(vectors) # scale the data - is this necessary/does it help or hinder? X = StandardScaler().fit_transform(X) # Compute DBSCAN params = self.params.multiple_lattice_search.cluster_analysis.dbscan db = DBSCAN(eps=params.eps, min_samples=params.min_samples).fit(X) core_samples = db.core_sample_indices_ # core_samples is a list of numpy.int64 objects core_samples = flex.int([int(i) for i in core_samples]) labels = flex.int(db.labels_.astype(np.int32)) # Number of clusters in labels, ignoring noise if present. n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0) logger.info("Estimated number of clusters: %d" % n_clusters_) return labels
[docs] def debug_plot_clusters(self, vectors, labels, min_cluster_size=1): assert len(vectors) == len(labels) from matplotlib import pyplot from mpl_toolkits.mplot3d import Axes3D # import dependency import numpy # Black removed and is used for noise instead. fig = pyplot.figure() ax = fig.add_subplot(111, projection="3d") # label clusters smaller than min_cluster_size as noise unique_labels = set(labels) for k in unique_labels: sel_k = labels == k if (sel_k).count(True) < min_cluster_size: labels = labels.set_selected(sel_k, -1) unique_labels = set(labels) n_clusters = len(unique_labels) - (1 if -1 in unique_labels else 0) colours = pyplot.cm.Spectral(numpy.linspace(0, 1, len(unique_labels))) # vectors = flex.vec3_double(vectors) ax.scatter([0], [0], [0], c="k", marker="+", s=50) for k, col in zip(unique_labels, colours): class_members = (labels == k).iselection() if k == -1: # or len(class_members) < min_cluster_size: # Black used for noise. col = "k" col = "0.25" # mid-grey markersize = 1 marker = "+" # continue else: markersize = 10 marker = "o" if not isinstance(col, basestring) and len(col) == 4: # darken the edges frac = 0.75 edgecolor = [col[0] * frac, col[1] * frac, col[2] * frac, col[3]] else: edgecolor = col vs = numpy.array([vectors[i] for i in class_members]) xs = vs[:, 0] ys = vs[:, 1] zs = vs[:, 2] # plot whole sphere for visual effect xs = numpy.concatenate((xs, -xs)) ys = numpy.concatenate((ys, -ys)) zs = numpy.concatenate((zs, -zs)) ax.scatter( xs, ys, zs, c=col, marker=marker, s=markersize, edgecolor=edgecolor ) pyplot.title("Estimated number of clusters: %d" % n_clusters) pyplot.show()
[docs] def find_peaks_clean(self): import omptbx # doesn't seem to be any benefit to using more than say 8 threads num_threads = min(8, omptbx.omp_get_num_procs(), self.params.nproc) omptbx.omp_set_num_threads(num_threads) d_min = self.params.fft3d.reciprocal_space_grid.d_min rlgrid = 2 / (d_min * self.gridding[0]) frame_number = self.reflections["xyzobs.px.value"].parts()[2] scan_range_min = max( int(math.floor(flex.min(frame_number))), self.imagesets[0].get_array_range()[0], ) # XXX what about multiple imagesets? scan_range_max = min( int(math.ceil(flex.max(frame_number))), self.imagesets[0].get_array_range()[1], ) # XXX what about multiple imagesets? scan_range = self.params.scan_range if not len(scan_range): scan_range = [[scan_range_min, scan_range_max]] scan = self.imagesets[0].get_scan() # XXX what about multiple imagesets? angle_ranges = [ [scan.get_angle_from_array_index(i, deg=False) for i in range_] for range_ in scan_range ] grid = flex.double(flex.grid(self.gridding), 0) sampling_volume_map( grid, flex.vec2_double(angle_ranges), self.imagesets[0].get_beam().get_s0(), self.imagesets[0].get_goniometer().get_rotation_axis(), rlgrid, d_min, self.params.b_iso, ) fft = fftpack.complex_to_complex_3d(self.gridding) grid_complex = flex.complex_double( reals=grid, imags=flex.double(grid.size(), 0) ) grid_transformed = fft.forward(grid_complex) grid_real = flex.pow2(flex.real(grid_transformed)) gamma = 1 peaks = flex.vec3_double() # n_peaks = 200 n_peaks = 100 # XXX how many do we need? dirty_beam = grid_real dirty_map = self.grid_real.deep_copy() import time t0 = time.time() peaks = clean_3d(dirty_beam, dirty_map, n_peaks, gamma=gamma) t1 = time.time() # print "clean_3d took %.2f s" %(t1-t0) reciprocal_lattice_points = self.reflections["rlp"].select( self.reflections_used_for_indexing ) peaks = self.optimise_peaks(peaks, reciprocal_lattice_points) peaks_frac = flex.vec3_double() for p in peaks: peaks_frac.append( ( p[0] / self.gridding[0], p[1] / self.gridding[1], p[2] / self.gridding[2], ) ) # print p, peaks_frac[-1] if self.params.debug: self.debug_write_ccp4_map(grid, "sampling_volume.map") self.debug_write_ccp4_map(grid_real, "sampling_volume_FFT.map") self.debug_write_ccp4_map(dirty_map, "clean.map") self.sites = peaks_frac # we don't really know the "volume"" of the peaks, but this method should # find the peaks in order of their intensity (strongest first) self.volumes = flex.double(range(len(self.sites), 0, -1)) return
[docs] def optimise_peaks(self, peaks, reciprocal_lattice_points): # optimise the peak position using a grid search around the starting peak position optimised_peaks = flex.vec3_double() n_points = 4 grid_step = 0.25 for peak in peaks: max_value = 1e-8 max_index = None for i in range(-n_points, n_points): i_coord = peak[0] + i * grid_step for j in range(-n_points, n_points): j_coord = peak[1] + j * grid_step for k in range(-n_points, n_points): k_coord = peak[2] + k * grid_step v = self.fft_cell.orthogonalize( ( i_coord / self.gridding[0], j_coord / self.gridding[1], k_coord / self.gridding[2], ) ) two_pi_S_dot_v = 2 * math.pi * reciprocal_lattice_points.dot(v) f = flex.sum(flex.cos(two_pi_S_dot_v)) if f > max_value: max_value = f max_index = (i_coord, j_coord, k_coord) optimised_peaks.append(max_index) return optimised_peaks
[docs]def sampling_volume_map( data, angle_range, beam_vector, rotation_axis, rl_grid_spacing, d_min, b_iso ): from dials.algorithms.indexing import sampling_volume_map return sampling_volume_map( data, angle_range, beam_vector, rotation_axis, rl_grid_spacing, d_min, b_iso )
[docs]def clean_3d(dirty_beam, dirty_map, n_peaks, gamma=1): from dials.algorithms.indexing import clean_3d as _clean_3d return _clean_3d(dirty_beam, dirty_map, n_peaks, gamma=gamma)