diff --git a/ica/python_dual_regress.py b/ica/python_dual_regress.py
index 73717af..e518712 100644
--- a/ica/python_dual_regress.py
+++ b/ica/python_dual_regress.py
@@ -135,7 +135,7 @@ def concat_regressors(a,b, outdir = None):
     return outf
         
 
-def sub_spatial_map(infile, design, mask, outdir, desnorm=True, mvt=None):
+def sub_spatial_map(infile, design, mask, outdir, desnorm=True, out_res=False, mvt=None):
     """ glm on ts data using stage1 txt file as model
     Parameters
     ----------
@@ -150,6 +150,8 @@ def sub_spatial_map(infile, design, mask, outdir, desnorm=True, mvt=None):
     desnorm : bool  (True, False, default = True)
         switch on normalisation of the design matrix
         columns to unit std. deviation
+    out_res : bool (True, False, default = False)
+        opion to output residuals image from glm
     mvt : file
         file containing movement regressors which will be concatenated
         to design (output from stage 1 glm)
@@ -159,18 +161,20 @@ def sub_spatial_map(infile, design, mask, outdir, desnorm=True, mvt=None):
         file of 4D timeseries components
     stage2_tsz : str
         file of z-transfomred 4D timesereis components
+    
 
     Notes
     -----
     fsl_glm -i <file> -d <outdir>/dr_stage1_${subid}.txt
     -o <outdir>/dr_stage2_$s --out_z=<outdir>/dr_stage2_${subid}_Z
-    --demean -m $OUTPUT/mask <desnorm>;
+    --demean -m $OUTPUT/mask <desnorm> --out_res=<outdir>/dr_stage2_${subid}_res;
     """
     
     subid = get_subid(infile)
     # define outfiles
     stage2_ts = os.path.join(outdir, 'dr_stage2_%s'%(subid))
     stage2_tsz = os.path.join(outdir,'dr_stage2_%s_Z'%(subid))
+    stage2_res = os.path.join(outdir,'dr_stage2_%s_res'%(subid))
     ext = get_fsl_outputtype()
     # add movment regressor to design if necessary
     if not mvt is None: 
@@ -185,6 +189,8 @@ def sub_spatial_map(infile, design, mask, outdir, desnorm=True, mvt=None):
     # Append des_norm flag to command if desnorm=True (defaults is True).
     if desnorm:
         cmd = ' '.join([cmd, '--des_norm'])
+    if out_res:
+        cmd = ' '.join([cmd, '--out_res=%s'%(stage2_res)])
     cout = CommandLine(cmd).run()
     if not cout.runtime.returncode == 0:
         print cmd
diff --git a/scripts/dualregress_userscript_example.py b/scripts/dualregress_userscript_example.py
index 4bebe4b..523a1c7 100644
--- a/scripts/dualregress_userscript_example.py
+++ b/scripts/dualregress_userscript_example.py
@@ -52,26 +52,32 @@
 if __name__ == '__main__':
 
     ###Set paths and specify ICA/template data
-    #########################################################
+    #######################################################################
     basedir = '/home/jagust/jelman/rsfmri_ica/data'
-
+    gica_dir = os.path.join(basedir, 'OldICA_IC0_ecat_2mm_6fwhm_125.gica')
+    # suffic of subject-specific ica directories created by melodic
+    # These will be appended to subj code later
+    sub_icadirname = '_4d_OldICA_IC0_ecat_2mm_6fwhm_125.ica' 
     ### get mask from groupica
-    mask = os.path.join(basedir, 'OldICA_IC40_ecat_7mm_125.gica', 
-                                        'groupmelodic.ica', 
-                                        'mask.nii.gz')
+    mask = os.path.join(gica_dir, 
+                        'groupmelodic.ica', 
+                        'mask.nii.gz')
 
     ### location for 4D template components (eg melodic_ICA, laird_4D)
-    template = os.path.join(basedir, 'OldICA_IC40_ecat_7mm_125.gica', 
-                                        'groupmelodic.ica', 
-                                        'melodic_IC.nii.gz')
-
+    template = os.path.join(gica_dir, 
+                            'groupmelodic.ica', 
+                            'melodic_IC.nii.gz')
+    ## Note: Set name of confound file name in the variable mvtfile located  
+    ## in dr stage 2 below.
+    #######################################################################
+    
     #Get number of vols in template file so that only ICs of interest
     #are merged, not those corresponding to nuisance regressors.
     cmd = ' '.join(['fslnvols', template])
     num_ics = int(commands.getoutput(cmd))
 
     ### output directory
-    outdir = os.path.join(basedir, 'OldICA_IC40_ecat_7mm_125.gica', 'dual_regress')
+    outdir = os.path.join(gica_dir, 'dual_regress')
     if os.path.isdir(outdir)==False:
         os.mkdir(outdir)      
     else:
@@ -83,7 +89,7 @@
         os.mkdir(outdir)
         print outdir, 'exists, moving to ', newdir
 
-    ### Specify input data filelist, otherwise searches basedir for data
+    ### Specify file listing input data files, otherwise searches basedir for data
     if len(sys.argv) ==2:   #If specified, load file as list of infiles
         args = sys.argv[1]
         print 'Using data specified in ', args
@@ -118,15 +124,15 @@
         ## If you want to add movement params & spike regressors to stage2 of model
         ## mvtfile = <confound file> 
         ## Must change input of mvt parameter from None in dr_stage2 below
-        sub_icadir = ''.join([subid, '_4d_OldICA_IC0_ecat_8mm_125.ica']) 
+        sub_icadir = ''.join([subid, sub_icadirname]) 
         mvtfile = os.path.join(basedir,
                             subid,
                             'func',
-                            sub_icadir, 'mc',
-                            'prefiltered_func_data_mcf.par') #Name of confound file
-
+                            'confound_regressors_6mm.txt') #Name of confound file
+                            
+        ## If you want residuals to be output, change out_res to True below
         stage2_ts, stage2_tsz = pydr.sub_spatial_map(tmpf, txtf, mask, outdir,
-                                     desnorm=True, mvt=mvtfile)
+                                     desnorm=True, out_res=True,mvt=mvtfile)
         
         subd.update({subid:stage2_ts})
 
diff --git a/tools/net_affiliation.py b/tools/net_affiliation.py
new file mode 100644
index 0000000..490cc03
--- /dev/null
+++ b/tools/net_affiliation.py
@@ -0,0 +1,100 @@
+import os, sys
+import numpy as np
+import nibabel as nib
+from glob import glob
+import general_utilities as gu
+
+def get_dims(data):
+    """
+    Given an array, returns dimensions necessary to broadcast to a 2D array.
+
+    Input:
+    data : numpy array 
+
+    Returns:
+    n_voxels : int
+            product of all axis shapes except the last
+    n_scans : int
+            equal to the shape of the last axis
+    """
+    n_voxels = np.prod(data.shape[:-1])
+    n_scans = data.shape[-1]
+    return n_voxels, n_scans
+
+def get_primary_net(dat_array2d):
+    """
+    Given a 2D numpy array, returns an array with the max value across
+    the last axis and a masked array with all values except the max
+
+    Input:
+    dat_array2d : numpy array
+
+    Returns:
+    prim_array : numpy array
+                array of max values across the last dimension
+    other_array : numpy masked array
+                masked array of all values from dat_array2d
+                except those contained in prim_array
+    """
+    assert len(dat_array2d.shape) == 2
+    prim_net_idx = dat_array2d.argmax(axis=1)
+    prim_net = dat_array2d.max(axis=1)
+    prim_net = np.reshape(prim_net, (dat_array2d.shape[0],1))
+    other_nets = np.ma.array(dat_array2d, mask=False)
+    other_nets.mask[np.arange(len(dat_array2d)),prim_net_idx] = True
+    return prim_net, other_nets
+
+def calculate_diff_map(dat_array):
+    """
+    Takes a numpy array and finds the mean difference between the max 
+    value and all others across time/scans. 
+
+    Input:
+    dat_array : numpy array (at least 2d)
+
+    Returns:
+    meandiff_array : numpy array
+                The mean difference between the primary network and
+                each of the other networks
+    sumdiff_array : numpy array
+                The sum of sifferences beterrn prmaru network and 
+                each of the other networks
+    splitdiff_array : numpy array
+                The difference between the primary network and the sum
+                of all other networks
+    """
+    dat_array = np.atleast_2d(dat_array)
+    n_voxels, n_nets = get_dims(dat_array)
+    dat_array2d = np.reshape(dat_array, (n_voxels, n_nets))
+    prim_net, other_nets = get_primary_net(dat_array2d) 
+    other_nets_mean = other_nets.mean(axis=1)
+    other_nets_mean = np.reshape(other_nets_mean, (other_nets_mean.shape[0],1))
+    nets_diff = prim_net - other_nets_mean
+    nets_diff_array = np.reshape(nets_diff, (dat_array.shape[:-1]))
+    return nets_diff_array
+
+def main(datafiles):
+    for subj_file in datafiles:
+        subid = gu.get_subid(subj_file, subjstr)
+        print 'Starting on subject %s'%(subid)
+        dat, aff = gu.load_nii(subj_file)
+        nets_dat = dat[:,:,:,net_idx]
+        nets_diff_array = calculate_diff_map(nets_dat)
+        nets_diff_outfile = os.path.join(outdir, ''.join(['nets_diff_',subid,'.nii.gz']))
+        gu.save_nii(nets_diff_outfile, nets_diff_array, aff)
+        
+        
+if __name__ == '__main__':
+        
+    #Path to 4D files in which networks are concatenated over time
+    datadir = '/home/jagust/rsfmri_ica/data/Allsubs_YoungICA_2mm_IC30.gica/dual_regress'
+    #Directory to output difference maps
+    outdir = '/home/jagust/rsfmri_ica/data/Allsubs_YoungICA_2mm_IC30.gica/difference_maps'
+    subjstr = 'subject[0-9]{5}'
+    dataglobstr = 'dr_stage2_subject*_Z.nii.gz'
+    datafiles = glob(os.path.join(datadir, dataglobstr))
+    #Indicies of networks to include. Start count at 0. 
+    net_idx = [0,1,2,3,4,6,7,8,9,12,14,15,24,29] 
+
+    main(datafiles)
+
diff --git a/tools/rapid_art.py b/tools/rapid_art.py
new file mode 100644
index 0000000..df70069
--- /dev/null
+++ b/tools/rapid_art.py
@@ -0,0 +1,208 @@
+import os, sys
+import numpy as np
+import nibabel as ni
+import nipype.algorithms.rapidart as rapidart
+from glob import glob
+import nipype.interfaces.fsl as fsl
+import nipype.algorithms.misc as misc
+import json
+import nipy
+import nipy.algorithms.diagnostics as diag
+from nipype.utils.filemanip import fname_presuffix
+import shutil
+import argparse
+
+def make_4d_nibabel(infiles,outdir=None):
+    shape = ni.load(infiles[0]).get_shape()
+    finalshape = tuple([x for x in shape]+[len(infiles)])
+    dat4d = np.empty(finalshape)
+    for val, f in enumerate(infiles):
+        tmpdat = ni.load(f).get_data()
+        tmpdat[np.isnan(tmpdat)] =  0
+        dat4d[:,:,:,val] = tmpdat
+    newimg = ni.Nifti1Image(dat4d, ni.load(infiles[0]).get_affine())
+    if outdir is None:
+        outf = fname_presuffix(infiles[0], prefix='data4d_')
+    else:
+        outf = fname_presuffix(infiles[0], prefix='data4d_',newpath = outdir)
+    newimg.to_filename(outf)
+    return outf
+
+        
+def make_4d(infiles):
+    merge = fsl.utils.Merge()
+    merge.inputs.in_files = infiles
+    merge.inputs.dimension = 't'
+    mrgout = merge.run()
+    return mrgout.outputs.merged_file
+    
+def outlier_vol_fromfile(mean, outliers, fname = None):
+    _, filename = os.path.split(mean)
+    vox = np.loadtxt(outliers).astype(int)
+    img = ni.load(mean)
+    newdat = np.zeros(img.get_shape()).flatten()
+    newdat[vox] = 100
+    newdat.shape = img.get_shape()
+    newimg = ni.Nifti1Image(newdat, img.get_affine())
+    if fname is None:
+        newfile = filename.replace('mean', 'outliers')
+    newfile = fname
+    newimg.to_filename(newfile)
+    print os.path.abspath(newfile)
+
+def gen_sig2noise_img(in4d, outdir):
+    startdir = os.getcwd()
+    os.chdir(outdir)
+    tsnr = misc.TSNR()
+    tsnr.inputs.in_file = in4d
+    tsnrout = tsnr.run()
+    ## there is no valid return code for tsnr 
+    #if tsnrout.runtime.returncode == 0:
+    #    return tsnrout.outputs.tsnr_file
+    #else:
+    #    print tsnrout.runtime.stderr
+    #    return None
+    os.chdir(startdir)
+    return tsnrout.outputs.tsnr_file
+
+def make_qa_dir(inroot, name='data_QA'):
+    qadir = os.path.join(inroot, name)
+    if os.path.isdir(qadir):
+        return qadir, True
+    else:
+        os.mkdir(qadir)
+        return qadir, False
+
+def run_artdetect(file4d, param_file, thresh = 4, param_source='SPM'):
+    startdir = os.getcwd()
+    pth, _ = os.path.split(param_file)
+    os.chdir(pth)
+    ad = rapidart.ArtifactDetect()
+    ad.inputs.realigned_files = file4d
+    ad.inputs.realignment_parameters = param_file
+    ad.inputs.parameter_source = param_source
+    ad.inputs.norm_threshold = 1
+    ad.inputs.mask_type = 'spm_global'
+    ad.inputs.use_differences = [True, False]
+    ad.inputs.zintensity_threshold = thresh
+    adout = ad.run()
+    os.chdir(startdir)
+    return adout
+
+def screen_data_dirnme(in4d, outdir):
+    """uses nipy diagnostic code to screen the data for
+    outlier values and saves results to three images
+    mean, std, pca, in same dir as original file(s)"""
+    img = nipy.load_image(in4d)
+    result = diag.screen(img)
+    # save mean, std, pca
+    pth, nme = os.path.split(in4d)
+    stripnme = nme.split('.')[0]
+    
+    pcafile = os.path.join(outdir,
+                           'QA-PCA_%s.nii.gz'%(nme))
+    meanfile = os.path.join(outdir,
+                            'QA-MEAN_%s.nii.gz'%(nme))
+    stdfile = os.path.join(outdir,
+                           'QA-STD_%s.nii.gz'%(nme))
+    nipy.save_image(result['mean'], meanfile)
+    nipy.save_image(result['std'], stdfile)
+    nipy.save_image(result['pca'], pcafile)
+    print 'saved: %s\n \t%s\n \t%s\n'%(pcafile, meanfile, stdfile)
+
+def save_qa_img_dirnme(in4d, outdir):
+    pth, nme = os.path.split(in4d)
+    img = nipy.load_image(in4d)
+    diag.plot_tsdiffs(diag.time_slice_diffs(img))
+    cleantime = time.asctime().replace(' ','-').replace(':', '_')
+    figfile = os.path.join(outdir, 'QA_%s_%s.png'%(nme, cleantime))
+    pylab.savefig(figfile)
+
+
+def main(infile, param_file, param_source, thresh, outdir=None):
+
+    
+    if len(infile) > 1:
+        # input is list, need to merge
+        if outdir is None:
+            outdir, _ = os.path.split(infile[0])
+        qadir, exists = make_qa_dir(outdir)
+        if exists:
+            print '%s exists, remove to re-run'%qadir
+            return None
+        merged = make_4d_nibabel(infile,outdir=qadir)
+    else:
+        merged = infile[0]
+        if outdir is None:
+            outdir, _ = os.path.split(merged)
+        qadir, exists = make_qa_dir(outdir)
+        if exists:
+            print '%s exists, remove to re-run'%qadir
+            return None
+    param_file = shutil.copyfile(param_file, qadir)
+    artout = run_artdetect(merged, param_file, thresh, param_source)
+
+    vox_outliers =  artout.outputs.outlier_files
+
+    statd = json.load(open(artout.outputs.statistic_files))
+    mot = statd[1]['motion_outliers']
+    intensity = statd[1]['intensity_outliers']
+    if mot > 0:
+        try:
+            tmp = np.loadtxt(artout.outputs.outlier_files)[:mot]
+        except:
+            tmp = np.loadtxt(artout.outputs.outlier_files)
+	tmp.tofile(os.path.join(qadir,'bad_movement_frames.txt'), 
+                   sep='\n')
+    np.array([mot,intensity]).tofile(os.path.join(qadir,
+                                                  'motion_intensity_outliers'),
+                                     sep = '\n')
+    print 'QA written to %s'%(qadir)
+        
+
+if __name__ == '__main__':
+
+    # create the parser
+    parser = argparse.ArgumentParser(
+        description='Artifact Detection on fMRI ')
+
+    # add the arguments
+    parser.add_argument(
+        'infiles',
+        type = str,
+        nargs = '+', # one or more items
+        help = 'file(s) 4d or multiple 3d of fMRI series')
+    parser.add_argument(
+        '-params',
+        dest = 'params',
+        help = 'Movement parameters txt file')
+    parser.add_argument(
+        '-params_source',
+        dest = 'params_source',
+        default = 'SPM',
+        help = 'Params are from SPM or FSL (defualt SPM), specify which')
+    parser.add_argument(
+        '-thresh',
+        type = int,
+        dest = 'thresh',
+        default = 4,
+        help = 'Intensity threshold (default 4)')
+    parser.add_argument(
+        '-outdir',        
+        dest = 'outdir',
+        help = """Optional output directory (QA_data folder will be created
+        in this base direcoty if sepcified,
+        otherwise in directory of infiles""")
+    
+    if len(sys.argv) ==1:
+        parser.print_help()
+    else:
+        args = parser.parse_args()
+        print args
+        main(args.infiles,
+             args.params,
+             args.params_source,
+             args.thresh,
+             args.outdir)
+        
+  
diff --git a/tools/run_image_qa.py b/tools/run_image_qa.py
index 81b5169..e657bc0 100755
--- a/tools/run_image_qa.py
+++ b/tools/run_image_qa.py
@@ -1,5 +1,4 @@
 import sys, os
-sys.path.append("/home/jagust/cindeem/CODE/PetProcessing/misc")
 import rapid_art
 import numpy as np
 
@@ -51,10 +50,10 @@ def CreateRegressors(funcdir, art_output, num_vols):
 
     Parameters
     -----------
-    outdir : str
-        directory to save regressor text file in
+    funcdir : str
+        subject functional directory containing QA folder
     art_output : str
-        path to directory of rapid_art.py output
+        name of rapid_art.py output containing bad volumes
     num_vols : int
         number of volumes in pre-processed functional file
     
@@ -72,21 +71,27 @@ def CreateRegressors(funcdir, art_output, num_vols):
     qa_file = os.path.join(funcdir,'data_QA',art_output)
     outliers = np.loadtxt(qa_file, dtype=int)
     outliers = np.atleast_1d(outliers)
-    if len(outliers) >= 1:
+    if len(outliers) > 1:
         exists = True
-        outlier_array = np.zeros(num_vols,dtype=float)
+        outlier_array = np.zeros((num_vols,len(outliers)),dtype=float)
         for i in range(len(outliers)):
-            outlier_array[outliers[i]]=1
+            outlier_array[outliers[i],i]=1
         outfile = os.path.join(funcdir, 'data_QA', 'outliers_for_fsl.txt')
         outlier_array.tofile(outfile)
         np.savetxt(outfile, outlier_array, fmt='%i', delimiter=u'\t')
         print 'Saved %s'%outfile
-        return exists, outlier_array
-    elif len(outliers)==0:
+    elif len(outliers) == 1:
+        exists = True
+        outlier_array = np.zeros((num_vols,len(outliers)),dtype=float)
+        outlier_array[outliers[0],0] = 1
+        outfile = os.path.join(funcdir, 'data_QA', 'outliers_for_fsl.txt')
+        outlier_array.tofile(outfile)
+        np.savetxt(outfile, outlier_array, fmt='%i', delimiter=u'\t')
+        print 'Saved %s'%outfile
+    else:
         outlier_array = np.array([])
         print 'No outliers, only mc parameters will be used'
-        return exists, outlier_array
-
+    return exists, outlier_array
 
 
 def CombineRegressors(mc_params, outlier_array, outdir, confound_outname):
@@ -148,6 +153,7 @@ def CombineRegressors(mc_params, outlier_array, outdir, confound_outname):
         outdir = funcdir
         confound_outname = 'confound_regressors_6mm.txt'
         ######################################################
+        
         #Run artdetect and create QA directory
         rapid_art.main(infiles, param_file, param_source, thresh, outdir)