code: from iconic_matcher import IconicMatcher
  #from realign4d import TimeSeries, realign4d, resample4d
  import transform
  
  from numpy.testing import Tester
  
  test = Tester().test
  bench = Tester().bench 
  
   #!/usr/bin/env python
  
  from nipy.testing import assert_equal, assert_almost_equal, assert_raises
  import numpy as np
  
  from nipy.neurospin.register.iconic_matcher import IconicMatcher
  
  
  class Image(object):
      """ 
      Empty object to easily create image objects independently from any I/O package.
      """
      def __init__(self, array, toworld=None, voxsize=[1, 1, 1]):
          self.array = array
          self.voxsize = np.asarray(voxsize)
          if toworld == None: 
              toworld = np.diag(np.concatenate((self.voxsize, [1])))
          self.toworld = toworld
  
  
  
  def make_data_uint8(dx=100, dy=100, dz=50):
      return (256*(np.random.rand(dx, dy, dz) - np.random.rand())).astype('uint8')
  
  def make_data_int16(dx=100, dy=100, dz=50):
      return (256*(np.random.rand(dx, dy, dz) - np.random.rand())).astype('int16')
  
  def make_data_float64(dx=100, dy=100, dz=50):
      return (256*(np.random.rand(dx, dy, dz) - np.random.rand())).astype('float64')
  
  def _test_clamping(I, thI=0.0, clI=256):
      IM = IconicMatcher(I.array, I.array, I.toworld, I.toworld, thI, thI, source_bins=clI, target_bins=clI)
      Ic = IM.source_clamped
      Ic2 = IM.target_clamped[1:I.array.shape[0]+1,1:I.array.shape[1]+1,1:I.array.shape[2]+1].squeeze()
      assert_equal(Ic, Ic2)
      dyn = Ic.max() + 1
      assert_equal(dyn, IM.joint_hist.shape[0])
      assert_equal(dyn, IM.joint_hist.shape[1])
      assert_equal(dyn, IM.source_hist.shape[0])
      assert_equal(dyn, IM.target_hist.shape[0])
  
  def test_clamping_uint8(): 
      I = Image(make_data_uint8())
      _test_clamping(I)
  
  def test_clamping_uint8_nonstd():
      I = Image(make_data_uint8())
      _test_clamping(I, 10, 165)
  
  def test_clamping_int16(): 
      I = Image(make_data_int16())
      _test_clamping(I)
  
  def test_clamping_int16_nonstd(): 
      I = Image(make_data_int16())
      _test_clamping(I, 10, 165)
  
  def test_clamping_float64(): 
      I = Image(make_data_float64())
      _test_clamping(I)
  
  def test_clamping_float64_nonstd():
      I = Image(make_data_float64())
      _test_clamping(I, 10, 165)
          
  def _test_similarity_measure(simi, val):
      I = Image(make_data_int16())
      J = Image(I.array.copy())
      IM = IconicMatcher(I.array, J.array, I.toworld, J.toworld)
      IM.set_field_of_view(subsampling=[2,1,3])
      IM.set_similarity(simi)
      assert_almost_equal(IM.eval(np.eye(4)), val)
  
  def test_correlation_coefficient():
      _test_similarity_measure('cc', 1.0) 
  
  def test_correlation_ratio():
      _test_similarity_measure('cr', 1.0) 
  
  def test_normalized_mutual_information():
      _test_similarity_measure('nmi', 1.0) 
  
  def test_explore(): 
      I = Image(make_data_int16())
      J = Image(make_data_int16())
      IM = IconicMatcher(I.array, J.array, I.toworld, J.toworld)
      T = np.eye(4)
      T[0:3,3] = np.random.rand(3)
      simi, params = IM.explore(ux=[-1,0,1],uy=[-1,0,1])
  
  def test_iconic():
      """ Test the iconic class.
      """
      I = Image(make_data_int16())
      J = Image(I.array.copy())
      IM = IconicMatcher(I.array, J.array, I.toworld, J.toworld)
      assert_raises(ValueError, IM.set_field_of_view, subsampling=[0,1,3])
  
  if __name__ == "__main__":
          import nose
          nose.run(argv=['', __file__])
   """
  Example of a script that crates a 'hierarchical roi' structure
  from the blob model of an image
  
  fixme : redo it
  Used mainly for debugging at the moment (before unittests are created)
  
  This example is based on a (simplistic) simulated image.
  # Author : Bertrand Thirion, 2008-2009
  """
  
  import numpy as np
  
  
  import nipy.neurospin.spatial_models.hroi as hroi
  import nipy.neurospin.utils.simul_2d_multisubject_fmri_dataset as simul
  import nipy.neurospin.graph.field as ff
  
  ##############################################################################
  # simulate the data
  dimx = 60
  dimy = 60
  pos = 2*np.array([[6,7],[10,10],[15,10]])
  ampli = np.array([3,4,4])
  
  dataset = simul.make_surrogate_array(nbsubj=1, dimx=dimx, dimy=dimy, pos=pos,
                                       ampli=ampli, width=10.0).squeeze()
  
  dataset = np.reshape(dataset, (dimx, dimy,1))
  ref_dim = (dimx,dimy,1)
  xyz = np.array(np.where(dataset)).T
  nbvox = np.size(xyz, 0)
  affine = np.eye(4)
  shape = (dimx,dimy,1)
  
  # create the field strcture that encodes image topology
  Fbeta = ff.Field(nbvox)
  Fbeta.from_3d_grid(xyz.astype(np.int), 18)
  beta = np.reshape(dataset,(nbvox,1))
  Fbeta.set_field(beta)
  nroi = hroi.NROI_from_field(Fbeta, affine, shape, xyz, th=2.0, smin = 10)
  if nroi != None:
      n1 = nroi.copy()
      n2 = nroi.reduce_to_leaves()
  
  td = n1.depth_from_leaves()
  a = np.argmax(td)
  lv = n1.rooted_subtree(a)
  u = nroi.cc()
  u = np.nonzero(u == u[0])[0]
  
  nroi.set_discrete_feature_from_index('activation',beta)
  nroi.discrete_to_roi_features('activation')
  
  
  label = -np.ones((dimx,dimy))
  for k in range(nroi.k):
      label[nroi.xyz[k][:,0],nroi.xyz[k][:,1]]=k
  import matplotlib.pylab as mp
  mp.figure()
  mp.imshow(label,interpolation='Nearest')
  mp.show()
  
apis: numpy.testing.Tester
apis: numpy.concatenate numpy.eye numpy.asarray numpy.random.rand
apis: numpy.eye numpy.ones numpy.reshape matplotlib.pylab.figure numpy.argmax numpy.size matplotlib.pylab.show numpy.array matplotlib.pylab.imshow numpy.where numpy.nonzero