# -----------------------------------------------------------------------------
# Example script for fitting one map in another without the graphical user
# interface.
#
#       chimera --nogui fitnogui.py
#
# It can also be run using the graphical Chimera interface using File/Open.
#
def fit_map_in_map(map1_path, map2_path,
                   initial_map1_transform = None,
                   map1_threshold = None,
                   ijk_step_size_min = 0.01,    # Grid index units
                   ijk_step_size_max = 0.5,     # Grid index units
                   max_steps = 100,
                   optimize_translation = True,
                   optimize_rotation = True):


  # Files have to have file suffix indicating volume format.
  from VolumeData import open_file
  data1 = open_file(map1_path)
  data2 = open_file(map2_path)

  from VolumeViewer import Data_Region, full_region, Rendering_Options
  ro = Rendering_Options()
  map1 = Data_Region(data1, full_region(data1.size), '', ro)
  map2 = Data_Region(data2, full_region(data2.size), '', ro)

  if initial_map1_transform:
    from PDBmatrices.matrices import chimera_xform
    xf = chimera_xform(initial_map1_transform)
    map1.surface_model().openState.globalXform(xf)
    
  use_threshold = (map1_threshold != None)
  
  from FitMap import map_points_and_weights, motion_to_maximum
  points, point_weights = map_points_and_weights(map1, use_threshold)

  if len(points) == 0:
    if use_threshold:
      print 'No grid points above map threshold.'
    else:
      print 'Map has no non-zero values.'
    return

  move_tf, stats = motion_to_maximum(points, point_weights, map2, max_steps,
                                     ijk_step_size_min, ijk_step_size_max,
                                     optimize_translation, optimize_rotation)

  if initial_map1_transform:
    from PDBmatrices.matrices import multiply_matrices
    move_tf = multiply_matrices(move_tf, initial_map1_transform)

  header = ('\nFit map %s in map %s using %d points\n'
            % (map1.full_name(), map2.full_name(), stats['points']) +
            '  correlation = %.4g, overlap = %.4g\n'
            % (stats['correlation'], stats['overlap']) +
            '  steps = %d, shift = %.3g, angle = %.3g degrees\n'
            % (stats['steps'], stats['shift'], stats['angle']))
  print header

  from FitMap.gui import transformation_description
  tfs = transformation_description(move_tf)
  print tfs

# -----------------------------------------------------------------------------
# Example call fitting map into itself.
#
p1 = p2 = '/usr/local/src/chimera-demos/volume/examples/groel.mrc'
fit_map_in_map(p1, p2)