| 1 | # -----------------------------------------------------------------------------
|
|---|
| 2 | # Example script for fitting one map in another without the graphical user
|
|---|
| 3 | # interface.
|
|---|
| 4 | #
|
|---|
| 5 | # chimera --nogui fitnogui.py
|
|---|
| 6 | #
|
|---|
| 7 | # It can also be run using the graphical Chimera interface using File/Open.
|
|---|
| 8 | #
|
|---|
| 9 | # The rotation and translation to perform the fit are output.
|
|---|
| 10 | #
|
|---|
| 11 | # Only a local optimization is done so the initial position must be close
|
|---|
| 12 | # to the correct fit.
|
|---|
| 13 | #
|
|---|
| 14 | def fit_map_in_map(map1_path, map2_path,
|
|---|
| 15 | initial_map1_transform = None,
|
|---|
| 16 | map1_threshold = None,
|
|---|
| 17 | ijk_step_size_min = 0.01, # Grid index units
|
|---|
| 18 | ijk_step_size_max = 0.5, # Grid index units
|
|---|
| 19 | max_steps = 100,
|
|---|
| 20 | optimize_translation = True,
|
|---|
| 21 | optimize_rotation = True):
|
|---|
| 22 |
|
|---|
| 23 |
|
|---|
| 24 | # Files have to have file suffix indicating volume format.
|
|---|
| 25 | from VolumeViewer import open_volume_file
|
|---|
| 26 | map1 = open_volume_file(map1_path)[0] # Assume files contain just one array
|
|---|
| 27 | map2 = open_volume_file(map2_path)[0]
|
|---|
| 28 |
|
|---|
| 29 | if initial_map1_transform:
|
|---|
| 30 | from Matrix import chimera_xform
|
|---|
| 31 | xf = chimera_xform(initial_map1_transform)
|
|---|
| 32 | map1.surface_model().openState.globalXform(xf)
|
|---|
| 33 |
|
|---|
| 34 | use_threshold = (map1_threshold != None)
|
|---|
| 35 |
|
|---|
| 36 | from FitMap.fitmap import map_points_and_weights, motion_to_maximum
|
|---|
| 37 | points, point_weights = map_points_and_weights(map1, use_threshold)
|
|---|
| 38 |
|
|---|
| 39 | if len(points) == 0:
|
|---|
| 40 | if use_threshold:
|
|---|
| 41 | print 'No grid points above map threshold.'
|
|---|
| 42 | else:
|
|---|
| 43 | print 'Map has no non-zero values.'
|
|---|
| 44 | return
|
|---|
| 45 |
|
|---|
| 46 | move_tf, stats = motion_to_maximum(points, point_weights, map2, max_steps,
|
|---|
| 47 | ijk_step_size_min, ijk_step_size_max,
|
|---|
| 48 | optimize_translation, optimize_rotation)
|
|---|
| 49 |
|
|---|
| 50 | import Matrix
|
|---|
| 51 | if initial_map1_transform:
|
|---|
| 52 | move_tf = Matrix.multiply_matrices(move_tf, initial_map1_transform)
|
|---|
| 53 |
|
|---|
| 54 | header = ('\nFit map %s in map %s using %d points\n'
|
|---|
| 55 | % (map1.name, map2.name, stats['points']) +
|
|---|
| 56 | ' correlation = %.4g, overlap = %.4g\n'
|
|---|
| 57 | % (stats['correlation'], stats['overlap']) +
|
|---|
| 58 | ' steps = %d, shift = %.3g, angle = %.3g degrees\n'
|
|---|
| 59 | % (stats['steps'], stats['shift'], stats['angle']))
|
|---|
| 60 | print header
|
|---|
| 61 |
|
|---|
| 62 | tfs = Matrix.transformation_description(move_tf)
|
|---|
| 63 | print tfs
|
|---|
| 64 |
|
|---|
| 65 | # -----------------------------------------------------------------------------
|
|---|
| 66 | # Example call fitting map into itself.
|
|---|
| 67 | #
|
|---|
| 68 | p1 = p2 = '/Users/goddard/Downloads/Chimera/EMDB/emd_1080.map'
|
|---|
| 69 | fit_map_in_map(p1, p2)
|
|---|