docs: def MRSvoxelStats(segfile, MRSfile=None, center=None, dim=None, subjID=None,
                    gareas=GAREAS,wareas=WAREAS,csfareas=CSFAREAS):
      """
      returns grey/white/CSF content within MRS voxel
  
      Parameters
      ----------
      segfile:  nifti file
          path to segmentation file with grey/white matter labels (freesurfer
          aseg file converted from mgz). 
  
      MRSfile: nifti file
          path to MRSfile of MRS voxel. provide either this or center + dim
  
      center : integer array
          [x,y,z] where x, y and z are the coordinates of the point of
          interest. Provide either MRSfile or center+dim 
      
      dim : float array
          dimensions of voxel in mm. Provide either MRSfile or center+dim
  
      subjID: string
          optional subject identifier. Defaults to nims scan number
  
      gareas, wareas, csfareas: arrays of integers
          arrays of freesurfer labels for gray, white, and csf areas
          respectively. Determines which areas are considered
          grey/white/csf. Some areas like brainstem were lumped with grey
          matter. Did not include hypointensities. 
     
      Returns
      -------
  
      raw number of voxels and proportions of grey, white, csf and
      non-grey-or-white-matter 
  
      Notes
      -----
      See http://miykael.github.io/nipype-beginner-s-guide/regionOfInterest.html
  
      """
      # subject ID
      if subjID==None:
          m=re.search('(\w+?)_*_',segfile)
          subjID=m.group(0)[:-1] 
  
      # get segmentation file
      aseg = nib.load(segfile)
      aseg_data = aseg.get_data().squeeze()
      aseg_aff = aseg.get_affine()
      segdir = os.path.dirname(segfile)
      segvoxdim=np.diagonal(aseg_aff)[:3]
  
      # get nifti file of MRS voxel if one is provided
      if MRSfile is not None:
          if center is not None or dim is not None:
              msg = 'provide EITHER MRSfile OR center and dim, not both!'
              raise ValueError(msg)
  
          mrs = nib.load(MRSfile)
          mrs_data = mrs.get_data().squeeze()
          mrs_aff = mrs.get_affine()        
          # This applies the concatenation of the transforms from mrs space to
          # the T1 space. [0,0,0] is the center of the MRS voxel:
          center = np.round(np.dot(np.dot(np.linalg.pinv(aseg_aff), mrs_aff),
                                   [0,0,0,1]))[:3].astype(int)
          
          dim=np.diagonal(mrs_aff)[:3]
      else: # no MRSfile
          if center==None or dim==None:
              msg = 'if no MRSfile is provided, provide center and '
              msg += 'dimensions of voxel'
              raise ValueError()    
  
      # calculate beginning corner of MRS voxel
      print('Creating mask with center: [%s, %s, %s]'%(center[0],
                                                      center[1],
                                                      str(center[2])))
  
      print('and dimensions: [%s, %s, %s] mm'%(dim[0], dim[1], dim[2]))
      
      # calculate roi mask with numpy
      # round up to nearest number of voxel units
      voxdim=np.zeros(3)
      voxdim[0]=np.ceil(np.abs(dim[0] / segvoxdim[0]))
      voxdim[1]=np.ceil(np.abs(dim[1] / segvoxdim[1]))
      voxdim[2]=np.ceil(np.abs(dim[2] / segvoxdim[2]))
  
      print ('MRS voxel dimensions in T1 voxel units: [%s, %s, %s]'%(voxdim[0],
                                                                    voxdim[1],
                                                                    voxdim[2]))
  
      # corners
      lcorner = [center[0] - voxdim[0] / 2,
                 center[1] - voxdim[1] / 2,
                 center[2] - voxdim[2] / 2]
  
      ucorner = [center[0] + voxdim[0] / 2,
                 center[1] + voxdim[1] / 2,
                 center[2] + voxdim[2] / 2]
      # create mask
      mdata = np.zeros(aseg_data.shape)
      for i in range(int(lcorner[0]),int(ucorner[0])):
          for j in range(int(lcorner[1]),int(ucorner[1])):
              for k in range(int(lcorner[2]),int(ucorner[2])):
                  mdata[i,j,k]=1
  
  
      # calculate grey/white/csf from freesurfer labels
      gdata = np.zeros(aseg_data.shape)
      wdata = np.zeros(aseg_data.shape)
      csfdata = np.zeros(aseg_data.shape)
  
      for data, areas in zip([gdata, wdata, csfdata], [gareas, wareas, csfareas]):
          for area in areas:
              data[np.where(aseg_data==area)] = 1
  
      # multiply voxel ROI with seg data
      gmasked= mdata * gdata
      wmasked= mdata * wdata
      csfmasked= mdata * csfdata
      
      # extract stats from a given segmentation
      total = np.sum(mdata)
      white = np.sum(wmasked)
      grey = np.sum(gmasked)
      csf = np.sum(csfmasked)
      nongmwm = total - grey - white
  
      # proportions
      pWhite = float(white) / total
      pGrey = float(grey) / total
      pCSF = float(csf) / total 
      pNongmwm = float(nongmwm) / total
  
      return (total, grey, white, csf,nongmwm, pGrey, pWhite,pCSF, pNongmwm)