content: # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
  # vi: set ft=python sts=4 ts=4 sw=4 et:
  """
  The rapidart module provides routines for artifact detection and region of
  interest analysis.
  
  These functions include:
  
    * ArtifactDetect: performs artifact detection on functional images
  
    * StimulusCorrelation: determines correlation between stimuli
      schedule and movement/intensity parameters
  
     Change directory to provide relative paths for doctests
     >>> import os
     >>> filepath = os.path.dirname( os.path.realpath( __file__ ) )
     >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data'))
     >>> os.chdir(datadir)
  """
  
  import os
  from copy import deepcopy
  from warnings import warn
  
  from nibabel import load, funcs, Nifti1Image
  import numpy as np
  from scipy import signal
  import scipy.io as sio
  from nipype.external import six
  
  from ..interfaces.base import (BaseInterface, traits, InputMultiPath,
                                      OutputMultiPath, TraitedSpec, File,
                                      BaseInterfaceInputSpec, isdefined)
  from ..utils.filemanip import filename_to_list, save_json, split_filename
  from ..utils.misc import find_indices
  
  from .. import logging, config
  iflogger = logging.getLogger('interface')
  
  
  def _get_affine_matrix(params, source):
      """Return affine matrix given a set of translation and rotation parameters
  
      params : np.array (upto 12 long) in native package format
      source : the package that generated the parameters
               supports SPM, AFNI, FSFAST, FSL, NIPY
      """
      if source == 'FSL':
          params = params[[3, 4, 5, 0, 1, 2]]
      elif source in ('AFNI', 'FSFAST'):
          params = params[np.asarray([4, 5, 3, 1, 2, 0]) + (len(params) > 6)]
          params[3:] = params[3:] * np.pi / 180.
      if source == 'NIPY':
          # nipy does not store typical euler angles, use nipy to convert
          from nipy.algorithms.registration import to_matrix44
          return to_matrix44(params)
      #process for FSL, SPM, AFNI and FSFAST
      rotfunc = lambda x: np.array([[np.cos(x), np.sin(x)],
                                    [-np.sin(x), np.cos(x)]])
      q = np.array([0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0])
      if len(params) < 12:
          params = np.hstack((params, q[len(params):]))
      params.shape = (len(params),)
      # Translation
      T = np.eye(4)
      T[0:3, -1] = params[0:3]
      # Rotation
      Rx = np.eye(4)
      Rx[1:3, 1:3] = rotfunc(params[3])
      Ry = np.eye(4)
      Ry[(0, 0, 2, 2), (0, 2, 0, 2)] = rotfunc(params[4]).ravel()
      Rz = np.eye(4)
      Rz[0:2, 0:2] = rotfunc(params[5])
      # Scaling
      S = np.eye(4)
      S[0:3, 0:3] = np.diag(params[6:9])
      # Shear
      Sh = np.eye(4)
      Sh[(0, 0, 1), (1, 2, 2)] = params[9:12]
      if source in ('AFNI', 'FSFAST'):
          return np.dot(T, np.dot(Ry, np.dot(Rx, np.dot(Rz, np.dot(S, Sh)))))
      return np.dot(T, np.dot(Rx, np.dot(Ry, np.dot(Rz, np.dot(S, Sh)))))
  
  
  def _calc_norm(mc, use_differences, source, brain_pts=None):
      """Calculates the maximum overall displacement of the midpoints
      of the faces of a cube due to translation and rotation.
  
      Parameters
      ----------
      mc : motion parameter estimates
          [3 translation, 3 rotation (radians)]
      use_differences : boolean
      brain_pts : [4 x n_points] of coordinates
  
      Returns
      -------
  
      norm : at each time point
      displacement : euclidean distance (mm) of displacement at each coordinate
  
      """
  
      if brain_pts is None:
          respos = np.diag([70, 70, 75])
          resneg = np.diag([-70, -110, -45])
          all_pts = np.vstack((np.hstack((respos, resneg)), np.ones((1, 6))))
          displacement = None
      else:
          all_pts = brain_pts
      n_pts = all_pts.size - all_pts.shape[1]
      newpos = np.zeros((mc.shape[0], n_pts))
      if brain_pts is not None:
          displacement = np.zeros((mc.shape[0], n_pts / 3))
      for i in range(mc.shape[0]):
          affine = _get_affine_matrix(mc[i, :], source)
          newpos[i, :] = np.dot(affine,
                                all_pts)[0:3, :].ravel()
          if brain_pts is not None:
              displacement[i, :] = \
                    np.sqrt(np.sum(np.power(np.reshape(newpos[i, :],
                                                       (3, all_pts.shape[1])) -
                                            all_pts[0:3, :],
                                            2),
                                   axis=0))
      # np.savez('displacement.npz', newpos=newpos, pts=all_pts)
      normdata = np.zeros(mc.shape[0])
      if use_differences:
          newpos = np.concatenate((np.zeros((1, n_pts)),
                                   np.diff(newpos, n=1, axis=0)), axis=0)
          for i in range(newpos.shape[0]):
              normdata[i] = \
               np.max(np.sqrt(np.sum(np.reshape(np.power(np.abs(newpos[i, :]), 2),
                                                (3, all_pts.shape[1])), axis=0)))
      else:
          newpos = np.abs(signal.detrend(newpos, axis=0, type='constant'))
          normdata = np.sqrt(np.mean(np.power(newpos, 2), axis=1))
      return normdata, displacement
  
  
  def _nanmean(a, axis=None):
      """Return the mean excluding items that are nan
  
      >>> a = [1, 2, np.nan]
      >>> _nanmean(a)
      1.5
  
      """
      if axis:
          return np.nansum(a, axis) / np.sum(1 - np.isnan(a), axis)
      else:
          return np.nansum(a) / np.sum(1 - np.isnan(a))
  
  
  class ArtifactDetectInputSpec(BaseInterfaceInputSpec):
      realigned_files = InputMultiPath(File(exists=True),
                                  desc="Names of realigned functional data files",
                                       mandatory=True)
      realignment_parameters = InputMultiPath(File(exists=True), mandatory=True,
                              desc=("Names of realignment parameters"
                                    "corresponding to the functional data files"))
      parameter_source = traits.Enum("SPM", "FSL", "AFNI", "NiPy", "FSFAST",
                                     desc="Source of movement parameters",
                                     mandatory=True)
      use_differences = traits.ListBool([True, False], minlen=2, maxlen=2,
                                        usedefault=True,
              desc=("Use differences between successive motion (first element)"
                    "and intensity paramter (second element) estimates in order"
                    "to determine outliers.  (default is [True, False])"))
      use_norm = traits.Bool(True, requires=['norm_threshold'],
                             desc=("Uses a composite of the motion parameters in "
                                   "order to determine outliers."),
                             usedefault=True)
      norm_threshold = traits.Float(desc=("Threshold to use to detect motion-rela"
                                          "ted outliers when composite motion is "
                                          "being used"), mandatory=True,
                                    xor=['rotation_threshold',
                                         'translation_threshold'])
      rotation_threshold = traits.Float(mandatory=True, xor=['norm_threshold'],
              desc=("Threshold (in radians) to use to detect rotation-related "
                    "outliers"))
      translation_threshold = traits.Float(mandatory=True, xor=['norm_threshold'],
              desc=("Threshold (in mm) to use to detect translation-related "
                   "outliers"))
      zintensity_threshold = traits.Float(mandatory=True,
              desc=("Intensity Z-threshold use to detection images that deviate "
                    "from the mean"))
      mask_type = traits.Enum('spm_global', 'file', 'thresh',
              desc=("Type of mask that should be used to mask the functional "
                    "data. *spm_global* uses an spm_global like calculation to "
                    "determine the brain mask. *file* specifies a brain mask "
                    "file (should be an image file consisting of 0s and 1s). "
                    "*thresh* specifies a threshold to use. By default all voxels"
                    "are used, unless one of these mask types are defined."),
                              mandatory=True)
      mask_file = File(exists=True,
                       desc="Mask file to be used if mask_type is 'file'.")
      mask_threshold = traits.Float(desc=("Mask threshold to be used if mask_type"
                                          " is 'thresh'."))
      intersect_mask = traits.Bool(True,
                                   desc=("Intersect the masks when computed from "
                                         "spm_global."))
      save_plot = traits.Bool(True, desc="save plots containing outliers",
                              usedefault=True)
      plot_type = traits.Enum('png', 'svg', 'eps', 'pdf',
                              desc="file type of the outlier plot",
                              usedefault=True)
      bound_by_brainmask = traits.Bool(False, desc=("use the brain mask to "
                                                   "determine bounding box"
                                                   "for composite norm (works"
                                                   "for SPM and Nipy - currently"
                                                   "inaccurate for FSL, AFNI"),
                                       usedefault=True)
      global_threshold = traits.Float(8.0, desc=("use this threshold when mask "
                                                 "type equal's spm_global"),
                                      usedefault=True)
  
  
  class ArtifactDetectOutputSpec(TraitedSpec):
      outlier_files = OutputMultiPath(File(exists=True),
              desc=("One file for each functional run containing a list of "
                    "0-based indices corresponding to outlier volumes"))
      intensity_files = OutputMultiPath(File(exists=True),
              desc=("One file for each functional run containing the global "
                    "intensity values determined from the brainmask"))
      norm_files = OutputMultiPath(File,
              desc=("One file for each functional run containing the composite "
                    "norm"))
      statistic_files = OutputMultiPath(File(exists=True),
              desc=("One file for each functional run containing information "
                    "about the different types of artifacts and if design info is"
                    " provided then details of stimulus correlated motion and a "
                    "listing or artifacts by event type."))
      plot_files = OutputMultiPath(File,
              desc=("One image file for each functional run containing the "
                    "detected outliers"))
      mask_files = OutputMultiPath(File,
              desc=("One image file for each functional run containing the mask"
                    "used for global signal calculation"))
      displacement_files = OutputMultiPath(File,
              desc=("One image file for each functional run containing the voxel"
                    "displacement timeseries"))
  
  
  class ArtifactDetect(BaseInterface):
      """Detects outliers in a functional imaging series
  
      Uses intensity and motion parameters to infer outliers. If `use_norm` is
      True, it computes the movement of the center of each face a cuboid centered
      around the head and returns the maximal movement across the centers.
  
  
      Examples
      --------
  
      >>> ad = ArtifactDetect()
      >>> ad.inputs.realigned_files = 'functional.nii'
      >>> ad.inputs.realignment_parameters = 'functional.par'
      >>> ad.inputs.parameter_source = 'FSL'
      >>> ad.inputs.norm_threshold = 1
      >>> ad.inputs.use_differences = [True, False]
      >>> ad.inputs.zintensity_threshold = 3
      >>> ad.run() # doctest: +SKIP
      """
  
      input_spec = ArtifactDetectInputSpec
      output_spec = ArtifactDetectOutputSpec
  
      def __init__(self, **inputs):
          super(ArtifactDetect, self).__init__(**inputs)
  
      def _get_output_filenames(self, motionfile, output_dir):
          """Generate output files based on motion filenames
  
          Parameters
          ----------
  
          motionfile: file/string
              Filename for motion parameter file
          output_dir: string
              output directory in which the files will be generated
          """
          if isinstance(motionfile, six.string_types):
              infile = motionfile
          elif isinstance(motionfile, list):
              infile = motionfile[0]
          else:
              raise Exception("Unknown type of file")
          _, filename, ext = split_filename(infile)
          artifactfile = os.path.join(output_dir, ''.join(('art.', filename,
                                                           '_outliers.txt')))
          intensityfile = os.path.join(output_dir, ''.join(('global_intensity.',
                                                            filename, '.txt')))
          statsfile = os.path.join(output_dir, ''.join(('stats.', filename,
                                                        '.txt')))
          normfile = os.path.join(output_dir, ''.join(('norm.', filename,
                                                       '.txt')))
          plotfile = os.path.join(output_dir, ''.join(('plot.', filename, '.',
                                                       self.inputs.plot_type)))
          displacementfile = os.path.join(output_dir, ''.join(('disp.',
                                                               filename, ext)))
          maskfile = os.path.join(output_dir, ''.join(('mask.', filename, ext)))
          return (artifactfile, intensityfile, statsfile, normfile, plotfile,
                  displacementfile, maskfile)
  
      def _list_outputs(self):
          outputs = self._outputs().get()
          outputs['outlier_files'] = []
          outputs['intensity_files'] = []
          outputs['statistic_files'] = []
          outputs['mask_files'] = []
          if isdefined(self.inputs.use_norm) and self.inputs.use_norm:
              outputs['norm_files'] = []
              if self.inputs.bound_by_brainmask:
                  outputs['displacement_files'] = []
          if isdefined(self.inputs.save_plot) and self.inputs.save_plot:
              outputs['plot_files'] = []
          for i, f in enumerate(filename_to_list(self.inputs.realigned_files)):
              (outlierfile, intensityfile, statsfile, normfile, plotfile,
               displacementfile, maskfile) = \
                                        self._get_output_filenames(f, os.getcwd())
              outputs['outlier_files'].insert(i, outlierfile)
              outputs['intensity_files'].insert(i, intensityfile)
              outputs['statistic_files'].insert(i, statsfile)
              outputs['mask_files'].insert(i, maskfile)
              if isdefined(self.inputs.use_norm) and self.inputs.use_norm:
                  outputs['norm_files'].insert(i, normfile)
                  if self.inputs.bound_by_brainmask:
                      outputs['displacement_files'].insert(i, displacementfile)
              if isdefined(self.inputs.save_plot) and self.inputs.save_plot:
                  outputs['plot_files'].insert(i, plotfile)
          return outputs
  
      def _plot_outliers_with_wave(self, wave, outliers, name):
          import matplotlib.pyplot as plt
          plt.plot(wave)
          plt.ylim([wave.min(), wave.max()])
          plt.xlim([0, len(wave) - 1])
          if len(outliers):
              plt.plot(np.tile(outliers[:, None], (1, 2)).T,
                       np.tile([wave.min(), wave.max()], (len(outliers), 1)).T,
                       'r')
          plt.xlabel('Scans - 0-based')
          plt.ylabel(name)
  
      def _detect_outliers_core(self, imgfile, motionfile, runidx, cwd=None):
          """
          Core routine for detecting outliers
          """
          if not cwd:
              cwd = os.getcwd()
  
          # read in functional image
          if isinstance(imgfile, six.string_types):
              nim = load(imgfile)
          elif isinstance(imgfile, list):
              if len(imgfile) == 1:
                  nim = load(imgfile[0])
              else:
                  images = [load(f) for f in imgfile]
                  nim = funcs.concat_images(images)
  
          # compute global intensity signal
          (x, y, z, timepoints) = nim.get_shape()
  
          data = nim.get_data()
          affine = nim.get_affine()
          g = np.zeros((timepoints, 1))
          masktype = self.inputs.mask_type
          if masktype == 'spm_global':  # spm_global like calculation
              iflogger.debug('art: using spm global')
              intersect_mask = self.inputs.intersect_mask
              if intersect_mask:
                  mask = np.ones((x, y, z), dtype=bool)
                  for t0 in range(timepoints):
                      vol = data[:, :, :, t0]
                      # Use an SPM like approach
                      mask_tmp = vol > \
                                 (_nanmean(vol) / self.inputs.global_threshold)
                      mask = mask * mask_tmp
                  for t0 in range(timepoints):
                      vol = data[:, :, :, t0]
                      g[t0] = _nanmean(vol[mask])
                  if len(find_indices(mask)) < (np.prod((x, y, z)) / 10):
                      intersect_mask = False
                      g = np.zeros((timepoints, 1))
              if not intersect_mask:
                  iflogger.info('not intersect_mask is True')
                  mask = np.zeros((x, y, z, timepoints))
                  for t0 in range(timepoints):
                      vol = data[:, :, :, t0]
                      mask_tmp = vol > \
                                   (_nanmean(vol) / self.inputs.global_threshold)
                      mask[:, :, :, t0] = mask_tmp
                      g[t0] = np.nansum(vol * mask_tmp)/np.nansum(mask_tmp)
          elif masktype == 'file':  # uses a mask image to determine intensity
              maskimg = load(self.inputs.mask_file)
              mask = maskimg.get_data()
              affine = maskimg.get_affine()
              mask = mask > 0.5
              for t0 in range(timepoints):
                  vol = data[:, :, :, t0]
                  g[t0] = _nanmean(vol[mask])
          elif masktype == 'thresh':  # uses a fixed signal threshold
              for t0 in range(timepoints):
                  vol = data[:, :, :, t0]
                  mask = vol > self.inputs.mask_threshold
                  g[t0] = _nanmean(vol[mask])
          else:
              mask = np.ones((x, y, z))
              g = _nanmean(data[mask > 0, :], 1)
  
          # compute normalized intensity values
          gz = signal.detrend(g, axis=0)  # detrend the signal
          if self.inputs.use_differences[1]:
              gz = np.concatenate((np.zeros((1, 1)), np.diff(gz, n=1, axis=0)),
                                  axis=0)
          gz = (gz - np.mean(gz)) / np.std(gz)  # normalize the detrended signal
          iidx = find_indices(abs(gz) > self.inputs.zintensity_threshold)
  
          # read in motion parameters
          mc_in = np.loadtxt(motionfile)
          mc = deepcopy(mc_in)
  
          (artifactfile, intensityfile, statsfile, normfile, plotfile,
           displacementfile, maskfile) = self._get_output_filenames(imgfile, cwd)
          mask_img = Nifti1Image(mask.astype(np.uint8), affine)
          mask_img.to_filename(maskfile)
  
          if self.inputs.use_norm:
              brain_pts = None
              if self.inputs.bound_by_brainmask:
                  voxel_coords = np.nonzero(mask)
                  coords = np.vstack((voxel_coords[0],
                                      np.vstack((voxel_coords[1],
                                                 voxel_coords[2])))).T
                  brain_pts = np.dot(affine,
                                     np.hstack((coords,
                                                np.ones((coords.shape[0], 1)))).T)
              # calculate the norm of the motion parameters
              normval, displacement = _calc_norm(mc,
                                                 self.inputs.use_differences[0],
                                                 self.inputs.parameter_source,
                                                 brain_pts=brain_pts)
              tidx = find_indices(normval > self.inputs.norm_threshold)
              ridx = find_indices(normval < 0)
              if displacement is not None:
                  dmap = np.zeros((x, y, z, timepoints), dtype=np.float)
                  for i in range(timepoints):
                      dmap[voxel_coords[0],
                           voxel_coords[1],
                           voxel_coords[2], i] = displacement[i, :]
                  dimg = Nifti1Image(dmap, affine)
                  dimg.to_filename(displacementfile)
          else:
              if self.inputs.use_differences[0]:
                  mc = np.concatenate((np.zeros((1, 6)),
                                       np.diff(mc_in, n=1, axis=0)),
                                      axis=0)
              traval = mc[:, 0:3]  # translation parameters (mm)
              rotval = mc[:, 3:6]  # rotation parameters (rad)
              tidx = find_indices(np.sum(abs(traval) >
                                         self.inputs.translation_threshold, 1)
                                  > 0)
              ridx = find_indices(np.sum(abs(rotval) >
                                         self.inputs.rotation_threshold, 1) > 0)
  
          outliers = np.unique(np.union1d(iidx, np.union1d(tidx, ridx)))
  
          # write output to outputfile
          np.savetxt(artifactfile, outliers, fmt='%d', delimiter=' ')
          np.savetxt(intensityfile, g, fmt='%.2f', delimiter=' ')
          if self.inputs.use_norm:
              np.savetxt(normfile, normval, fmt='%.4f', delimiter=' ')
  
          if isdefined(self.inputs.save_plot) and self.inputs.save_plot:
              import matplotlib
              matplotlib.use(config.get("execution", "matplotlib_backend"))
              import matplotlib.pyplot as plt
              fig = plt.figure()
              if isdefined(self.inputs.use_norm) and self.inputs.use_norm:
                  plt.subplot(211)
              else:
                  plt.subplot(311)
              self._plot_outliers_with_wave(gz, iidx, 'Intensity')
              if isdefined(self.inputs.use_norm) and self.inputs.use_norm:
                  plt.subplot(212)
                  self._plot_outliers_with_wave(normval, np.union1d(tidx, ridx),
                                                'Norm (mm)')
              else:
                  diff = ''
                  if self.inputs.use_differences[0]:
                      diff = 'diff'
                  plt.subplot(312)
                  self._plot_outliers_with_wave(traval, tidx,
                                                'Translation (mm)' + diff)
                  plt.subplot(313)
                  self._plot_outliers_with_wave(rotval, ridx,
                                                'Rotation (rad)' + diff)
              plt.savefig(plotfile)
              plt.close(fig)
  
          motion_outliers = np.union1d(tidx, ridx)
          stats = [{'motion_file': motionfile,
                    'functional_file': imgfile},
                   {'common_outliers': len(np.intersect1d(iidx, motion_outliers)),
                    'intensity_outliers': len(np.setdiff1d(iidx,
                                                           motion_outliers)),
                    'motion_outliers': len(np.setdiff1d(motion_outliers, iidx)),
                    },
           {'motion': [{'using differences': self.inputs.use_differences[0]},
                        {'mean': np.mean(mc_in, axis=0).tolist(),
                         'min': np.min(mc_in, axis=0).tolist(),
                         'max': np.max(mc_in, axis=0).tolist(),
                         'std': np.std(mc_in, axis=0).tolist()},
                        ]},
           {'intensity': [{'using differences': self.inputs.use_differences[1]},
                           {'mean': np.mean(gz, axis=0).tolist(),
                            'min': np.min(gz, axis=0).tolist(),
                            'max': np.max(gz, axis=0).tolist(),
                            'std': np.std(gz, axis=0).tolist()},
                           ]},
                   ]
          if self.inputs.use_norm:
              stats.insert(3, {'motion_norm':
                                   {'mean': np.mean(normval, axis=0).tolist(),
                                    'min': np.min(normval, axis=0).tolist(),
                                    'max': np.max(normval, axis=0).tolist(),
                                    'std': np.std(normval, axis=0).tolist(),
                                    }})
          save_json(statsfile, stats)
  
      def _run_interface(self, runtime):
          """Execute this module.
          """
          funcfilelist = filename_to_list(self.inputs.realigned_files)
          motparamlist = filename_to_list(self.inputs.realignment_parameters)
          for i, imgf in enumerate(funcfilelist):
              self._detect_outliers_core(imgf, motparamlist[i], i,
                                         cwd=os.getcwd())
          return runtime
  
  
  class StimCorrInputSpec(BaseInterfaceInputSpec):
      realignment_parameters = InputMultiPath(File(exists=True), mandatory=True,
          desc=('Names of realignment parameters corresponding to the functional '
                'data files'))
      intensity_values = InputMultiPath(File(exists=True), mandatory=True,
                desc='Name of file containing intensity values')
      spm_mat_file = File(exists=True, mandatory=True,
                          desc='SPM mat file (use pre-estimate SPM.mat file)')
      concatenated_design = traits.Bool(mandatory=True,
                desc='state if the design matrix contains concatenated sessions')
  
  
  class StimCorrOutputSpec(TraitedSpec):
      stimcorr_files = OutputMultiPath(File(exists=True),
                       desc='List of files containing correlation values')
  
  
  class StimulusCorrelation(BaseInterface):
      """Determines if stimuli are correlated with motion or intensity
      parameters.
  
      Currently this class supports an SPM generated design matrix and requires
      intensity parameters. This implies that one must run
      :ref:`ArtifactDetect <nipype.algorithms.rapidart.ArtifactDetect>`
      and :ref:`Level1Design <nipype.interfaces.spm.model.Level1Design>` prior to running this or
      provide an SPM.mat file and intensity parameters through some other means.
  
      Examples
      --------
  
      >>> sc = StimulusCorrelation()
      >>> sc.inputs.realignment_parameters = 'functional.par'
      >>> sc.inputs.intensity_values = 'functional.rms'
      >>> sc.inputs.spm_mat_file = 'SPM.mat'
      >>> sc.inputs.concatenated_design = False
      >>> sc.run() # doctest: +SKIP
  
      """
  
      input_spec = StimCorrInputSpec
      output_spec = StimCorrOutputSpec
  
      def _get_output_filenames(self, motionfile, output_dir):
          """Generate output files based on motion filenames
  
          Parameters
          ----------
          motionfile: file/string
              Filename for motion parameter file
          output_dir: string
              output directory in which the files will be generated
          """
          (_, filename) = os.path.split(motionfile)
          (filename, _) = os.path.splitext(filename)
          corrfile = os.path.join(output_dir, ''.join(('qa.', filename,
                                                       '_stimcorr.txt')))
          return corrfile
  
      def _stimcorr_core(self, motionfile, intensityfile, designmatrix, cwd=None):
          """
          Core routine for determining stimulus correlation
  
          """
          if not cwd:
              cwd = os.getcwd()
          # read in motion parameters
          mc_in = np.loadtxt(motionfile)
          g_in = np.loadtxt(intensityfile)
          g_in.shape = g_in.shape[0], 1
          dcol = designmatrix.shape[1]
          mccol = mc_in.shape[1]
          concat_matrix = np.hstack((np.hstack((designmatrix, mc_in)), g_in))
          cm = np.corrcoef(concat_matrix, rowvar=0)
          corrfile = self._get_output_filenames(motionfile, cwd)
          # write output to outputfile
          file = open(corrfile, 'w')
          file.write("Stats for:\n")
          file.write("Stimulus correlated motion:\n%s\n" % motionfile)
          for i in range(dcol):
              file.write("SCM.%d:" % i)
              for v in cm[i, dcol + np.arange(mccol)]:
                  file.write(" %.2f" % v)
              file.write('\n')
          file.write("Stimulus correlated intensity:\n%s\n" % intensityfile)
          for i in range(dcol):
              file.write("SCI.%d: %.2f\n" % (i, cm[i, -1]))
          file.close()
  
      def _get_spm_submatrix(self, spmmat, sessidx, rows=None):
          """
          Parameters
          ----------
          spmmat: scipy matlab object
              full SPM.mat file loaded into a scipy object
          sessidx: int
              index to session that needs to be extracted.
          """
          designmatrix = spmmat['SPM'][0][0].xX[0][0].X
          U = spmmat['SPM'][0][0].Sess[0][sessidx].U[0]
          if rows is None:
              rows = spmmat['SPM'][0][0].Sess[0][sessidx].row[0] - 1
          cols = spmmat['SPM'][0][0].Sess[0][sessidx].col[0][range(len(U))] - 1
          outmatrix = designmatrix.take(rows.tolist(), axis=0).take(cols.tolist(),
                                                                    axis=1)
          return outmatrix
  
      def _run_interface(self, runtime):
          """Execute this module.
          """
          motparamlist = self.inputs.realignment_parameters
          intensityfiles = self.inputs.intensity_values
          spmmat = sio.loadmat(self.inputs.spm_mat_file, struct_as_record=False)
          nrows = []
          for i in range(len(motparamlist)):
              sessidx = i
              rows = None
              if self.inputs.concatenated_design:
                  sessidx = 0
                  mc_in = np.loadtxt(motparamlist[i])
                  rows = np.sum(nrows) + np.arange(mc_in.shape[0])
                  nrows.append(mc_in.shape[0])
              matrix = self._get_spm_submatrix(spmmat, sessidx, rows)
              self._stimcorr_core(motparamlist[i], intensityfiles[i],
                                  matrix, os.getcwd())
          return runtime
  
      def _list_outputs(self):
          outputs = self._outputs().get()
          files = []
          for i, f in enumerate(self.inputs.realignment_parameters):
              files.insert(i, self._get_output_filenames(f, os.getcwd()))
          if files:
              outputs['stimcorr_files'] = files
          return outputs
  