class SetResolution:
    def __init__(self, filesToResample, resolution):
        """During initialization make sure all files are resampled
           at resolution we'd like to use for each pipeline stage
        """
        self.p = Pipeline()
        
        for FH in filesToResample:
            dirForOutput = self.getOutputDirectory(FH)
            currentRes = volumeFromFile(FH.getLastBasevol()).separations
            if not abs(abs(currentRes[0]) - abs(resolution)) < 0.01:
                crop = ma.autocrop(resolution, FH, defaultDir=dirForOutput)
                self.p.addStage(crop)
                mask = FH.getMask()
                if mask:
                    #Need to resample the mask as well.
                    cropMask = ma.mincresampleMask(FH,
                                                   FH,
                                                   outputLocation=FH,
                                                   likeFile=FH)
                    self.p.addStage(cropMask)
        
    def getOutputDirectory(self, FH):
        """Sets output directory based on whether or not we have a full
        RegistrationPipeFH class or we are just using RegistrationFHBase"""
        if isinstance(FH, rfh.RegistrationPipeFH):
            outputDir = "resampled"
        else:
            outputDir = FH.basedir
        return outputDir

def resampleToCommon(xfm, FH, statsGroup, statsKernels, nlinFH):
    blurs = []
    if isinstance(statsKernels, list):
        blurs = statsKernels
    elif isinstance(statsKernels, str):
        for i in statsKernels.split(","):
            blurs.append(float(i))
    else:
        print("Improper type of blurring kernels specified for stats calculation: " + str(statsKernels))
        sys.exit()
    pipeline = Pipeline()
    outputDirectory = FH.statsDir
    filesToResample = []
    for b in blurs:
        filesToResample.append(statsGroup.relativeJacobians[b])
        if statsGroup.absoluteJacobians:
            filesToResample.append(statsGroup.absoluteJacobians[b])
    for f in filesToResample:
        outputBase = removeBaseAndExtension(f).split(".mnc")[0]
        outputFile = createBaseName(outputDirectory, outputBase + "_common" + ".mnc")
        logFile = fh.logFromFile(FH.logDir, outputFile)
        targetAndLike=nlinFH.getLastBasevol()
        res = ma.mincresample(f, 
                              targetAndLike,
                              likeFile=targetAndLike,
                              transform=xfm,
                              output=outputFile,
                              logFile=logFile,
                              argArray=["-sinc"]) 
        pipeline.addStage(res)
    
    return pipeline

def maskFiles(FH, isAtlas, numAtlases=1):
    """ Assume that if there is more than one atlas, multiple
        masks were generated and we need to perform a voxel_vote. 
        Otherwise, assume we are using inputLabels from crossing with
        only one atlas. 
    """
    #MF TODO: Make this more general to handle pairwise option. 
    p = Pipeline()
    if not isAtlas:
        if numAtlases > 1:
            voxel = voxelVote(FH, False, True)
            p.addStage(voxel)
            mincMathInput = voxel.outputFiles[0]  
        else:
            mincMathInput = FH.returnLabels(True)[0]
        FH.setMask(mincMathInput)
    else:
        mincMathInput = FH.getMask()
    mincMathOutput = fh.createBaseName(FH.resampledDir, FH.basename)
    mincMathOutput += "_masked.mnc"   
    logFile = fh.logFromFile(FH.logDir, mincMathOutput)
    cmd = ["mincmath"] + ["-clobber"] + ["-mult"]
    cmd += [InputFile(mincMathInput)] + [InputFile(FH.getLastBasevol())] 
    cmd += [OutputFile(mincMathOutput)]
    mincMath = CmdStage(cmd)
    mincMath.setLogFile(LogFile(logFile))
    p.addStage(mincMath)
    FH.setLastBasevol(mincMathOutput)
    return(p)

class LabelAndFileResampling:
    def __init__(self, 
                 inputPipeFH,
                 templatePipeFH, 
                 name="initial", 
                 createMask=False):  
        self.p = Pipeline()
        self.name = name
        
        if createMask:
            resampleDefault = "tmp"
            labelsDefault = "tmp"
        else:
            resampleDefault = "resampled"
            labelsDefault = "labels"
              
        # Resample all inputLabels 
        inputLabelArray = templatePipeFH.returnLabels(True)
        if len(inputLabelArray) > 0:
            """ for the initial registration, resulting labels should be added
                to inputLabels array for subsequent pairwise registration
                otherwise labels should be added to labels array for voting """
            if self.name == "initial":
                addOutputToInputLabels = True
            else:
                addOutputToInputLabels = False
            for i in range(len(inputLabelArray)):
                """Note: templatePipeFH and inputPipeFH have the reverse order
                   from how they are passed into this function. This is intentional
                   because the mincresample classes use the first argument as the 
                   one from which to get the file to be resampled. Here, either the 
                   mask or labels to be resampled come from the template."""
                if createMask:
                    resampleStage = ma.mincresampleMask(templatePipeFH,
                                                        inputPipeFH,
                                                        defaultDir=labelsDefault,
                                                        likeFile=inputPipeFH,
                                                        argArray=["-invert"],
                                                        outputLocation=inputPipeFH,
                                                        labelIndex=i,
                                                        setInputLabels=addOutputToInputLabels)
                else:
                    resampleStage = ma.mincresampleLabels(templatePipeFH,
                                                          inputPipeFH, 
                                                          defaultDir=labelsDefault,
                                                          likeFile=inputPipeFH,
                                                          argArray=["-invert"],
                                                          outputLocation=inputPipeFH,
                                                          labelIndex=i,
                                                          setInputLabels=addOutputToInputLabels)
                self.p.addStage(resampleStage)
            # resample files
            resampleStage = ma.mincresample(templatePipeFH,
                                            inputPipeFH,
                                            defaultDir=resampleDefault,
                                            likeFile=inputPipeFH,
                                            argArray=["-invert"],
                                            outputLocation=inputPipeFH)
            self.p.addStage(resampleStage)

def getXfms(nlinFH, subjects, space, mbmDir, time=None):

    """For each file in the build-model registration (associated with the specified
       time point), do the following:
       
       1. Find the to-native.xfm for that file. 
       2. Find the matching subject at the specified time point
       3. Set this xfm to be the last xfm from nlin average to subject from step #2. 
       4. Find the -from-native.xfm file.
       5. Set this xfm to be the last xfm from subject to nlin.
       
       Note: assume that the names in processedDir match beginning file 
             names for each subject
             We are also assuming subjects is either a dictionary or a list. 
    """
    
    """First handle subjects if dictionary or list"""
    if isinstance(subjects, list):
        inputs = subjects
    elif isinstance(subjects, dict):
        inputs = []
        for s in subjects:
            inputs.append(subjects[s][time])
    else:
        logger.error("getXfms only takes a dictionary or list of subjects. Incorrect type has been passed. Exiting...")
        sys.exit()
    
    pipeline = Pipeline()
    baseNames = walk(mbmDir).next()[1]
    for b in baseNames:
        if space == "lsq6":
            xfmToNative = abspath(mbmDir + "/" + b + "/transforms/" + b + "-final-to_lsq6.xfm")
        elif space == "lsq12":
            xfmToNative = abspath(mbmDir + "/" + b + "/transforms/" + b + "-final-nlin.xfm")
            xfmFromNative = abspath(mbmDir + "/" + b + "/transforms/" + b + "_inv_nonlinear.xfm")
        elif space == "native":
            xfmToNative = abspath(mbmDir + "/" + b + "/transforms/" + b + "-to-native.xfm")
            xfmFromNative = abspath(mbmDir + "/" + b + "/transforms/" + b + "-from-native.xfm")
        else:
            logger.error("getXfms can only retrieve transforms to and from native, lsq6 or lsq12 space. Invalid parameter has been passed.")
            sys.exit()
        for inputFH in inputs:
            if fnmatch.fnmatch(inputFH.getLastBasevol(), "*" + b + "*"):
                if space=="lsq6":
                    ix = ma.xfmInvert(xfmToNative, inputFH)
                    pipeline.addStage(ix)
                    xfmFromNative = ix.outputFiles[0]
                nlinFH.setLastXfm(inputFH, xfmToNative)
                inputFH.setLastXfm(nlinFH, xfmFromNative)
    return pipeline

    def __init__(self,
                 inputFH,
                 targetFH, 
                 blurs=[0.3, 0.2, 0.15], 
                 step=[1,0.5,0.333333333333333],
                 gradient=[False,True,False],
                 simplex=[3,1.5,1],
                 defaultDir="tmp"):                                      

        # TO DO: Might want to take this out and pass in # of generations, since
        # checking happens there. 
        if len(blurs) == len(step) == len(simplex):
            # do nothing - all lengths are the same and we're therefore happy
            pass
        else:
            logger.error("The same number of entries are required for blurs, step, and simplex in LSQ12")
            sys.exit()
                
        self.p = Pipeline()
        self.inputFH = inputFH
        self.targetFH = targetFH
        self.blurs = blurs
        self.step = step
        self.blurs = blurs
        self.gradient = gradient
        self.simplex = simplex
        self.defaultDir = defaultDir
            
        self.blurFiles()
        self.buildPipeline()

    def __init__(self, 
                 inputArray, 
                 targetFH, 
                 nlinOutputDir, 
                 avgPrefix, 
                 nlin_protocol,
                 resolution=None):
        self.p = Pipeline()
        """Initial inputs should be an array of fileHandlers with lastBasevol in lsq12 space"""
        self.inputs = inputArray
        """Initial target should be the file handler for the lsq12 average"""
        self.target = targetFH 
        """Output directory should be _nlin """
        self.nlinDir = nlinOutputDir
        """Prefix to pre-pend to averages at each generation"""
        self.avgPrefix = avgPrefix
        """Empty array that we will fill with averages as we create them"""
        self.nlinAverages = [] 
        """Create the blurring resolution from the file resolution"""
        self.fileRes = resolution
        # hack:
        self.generations = 0

        if (nlin_protocol==None and resolution == None):
            print("\nError: NLIN module was initialized without a protocol, and without a resolution for the registrations to be run at. Please specify one of the two. Exiting\n") 
            sys.exit()
        if (nlin_protocol and resolution):
            # we should have the nlin_protocol be able to overwrite the given resolution:
            self.fileRes = None
        
        # Create new nlin group for each input prior to registration
        for i in range(len(self.inputs)):
            self.inputs[i].newGroup(groupName="nlin")

 def __init__(self, 
               targetOutputDir, #Output directory for files related to initial target (often _lsq12)
               inputFiles, 
               nlinDir, 
               avgPrefix, #Prefix for nlin-1.mnc, ... nlin-k.mnc 
               createAvg=True, #True=call mincAvg, False=targetAvg already exists
               targetAvg=None, #Optional path to initial target - passing name does not guarantee existence
               targetMask=None, #Optional path to mask for initial target
               nlin_protocol=None,
               reg_method=None):
     self.p = Pipeline()
     self.targetOutputDir = targetOutputDir
     self.inputFiles = inputFiles
     self.nlinDir = nlinDir
     self.avgPrefix = avgPrefix
     self.createAvg = createAvg
     self.targetAvg = targetAvg
     self.targetMask = targetMask
     self.nlin_protocol = nlin_protocol
     self.reg_method = reg_method
     
     # setup initialTarget (if needed) and initialize non-linear module
     self.setupTarget()
     self.initNlinModule()
     
     #iterate through non-linear registration and setup averages
     self.nlinModule.iterate()
     self.p.addPipeline(self.nlinModule.p)
     self.nlinAverages = self.nlinModule.nlinAverages
     self.nlinParams = self.nlinModule.nlinParams

def MAGeTRegister(inputFH, 
                  templateFH, 
                  regMethod,
                  name="initial", 
                  createMask=False,
                  lsq12_protocol=None,
                  nlin_protocol=None):
    
    p = Pipeline()
    if createMask:
        defaultDir="tmp"
    else:
        defaultDir="transforms"
    if regMethod == "minctracc":
        sp = HierarchicalMinctracc(inputFH, 
                                   templateFH,
                                   lsq12_protocol=lsq12_protocol,
                                   nlin_protocol=nlin_protocol,
                                   defaultDir=defaultDir)
        p.addPipeline(sp.p)
    elif regMethod == "mincANTS":
        register = LSQ12ANTSNlin(inputFH, 
                                 templateFH, 
                                 lsq12_protocol=lsq12_protocol,
                                 nlin_protocol=nlin_protocol,
                                 defaultDir=defaultDir)
        p.addPipeline(register.p)
        
    rp = LabelAndFileResampling(inputFH, templateFH, name=name, createMask=createMask)
    p.addPipeline(rp.p)
    
    return(p)

 def __init__(self, 
              inputPipeFH, 
              templatePipeFH,
              blurs=[1, 0.5, 0.3]):
     
     self.p = Pipeline()
     self.inputPipeFH = inputPipeFH
     self.templatePipeFH = templatePipeFH
     
     self.blurFiles(blurs)

 def __init__(self, inputs, dirs, options, avgPrefix=None, initModel=None):
     self.inputs = inputs
     self.dirs = dirs
     self.options = options
     self.avgPrefix = avgPrefix
     self.initModel = initModel
     self.nlinFH = None
     
     self.p = Pipeline()
     
     self.buildPipeline()

    def __init__(self,
                 inFile,
                 targetFile,
                 nameForStage=None,
                 **kwargs):
        self.p = Pipeline()
        self.outputFiles = [] # this will contain the outputFiles from the mincresample of the main MINC file
        self.outputFilesMask = [] # this will contain the outputFiles from the mincresample of the mask belonging to the main MINC file
        

        # the first step is to simply run the mincresample command:
        fileRS = mincresample(inFile,
                              targetFile,
                              **kwargs)
        if(nameForStage):
            fileRS.name = nameForStage
        self.p.addStage(fileRS)
        self.outputFiles = fileRS.outputFiles
        
        # initialize the array of outputs for the mask in case there is none to be resampled
        self.outputFilesMask = [None] * len(self.outputFiles)
        
        # next up, is this a file handler, and if so is there a mask that needs to be resampled?
        if(isFileHandler(inFile)):
            if(inFile.getMask()):
                # there is a mask associated with this file, should be updated
                # we have to watch out in terms of interpolation arguments, if 
                # the original resample command asked for "-sinc" or "-tricubic"
                # for instance, we should remove that argument for the mask resampling
                # these options would reside in the argArray...
                maskArgs = copy.deepcopy(kwargs)
                if(maskArgs["argArray"]):
                    argList = maskArgs["argArray"]
                    for i in range(len(argList)):
                        if(re.match("-sinc", argList[i]) or
                           re.match("-trilinear", argList[i]) or
                           re.match("-tricubic", argList[i]) ):
                            del argList[i]
                    maskArgs["argArray"] = argList                   
                
                # if the output file for the mincresample command was already
                # specified, add "_mask.mnc" to it
                if(maskArgs["output"]):
                    maskArgs["output"] = re.sub(".mnc", "_mask.mnc", maskArgs["output"])
                    
                maskRS = mincresampleMask(inFile,
                                          targetFile,
                                          **maskArgs)
                if(nameForStage):
                    maskRS.name = nameForStage + "--mask--"
                self.p.addStage(maskRS)   
                self.outputFilesMask = maskRS.outputFiles

 def __init__(self, inputArray, targetFH, nlinOutputDir, nlin_protocol=None):
     self.p = Pipeline()
     """Initial inputs should be an array of fileHandlers with lastBasevol in lsq12 space"""
     self.inputs = inputArray
     """Initial target should be the file handler for the lsq12 average"""
     self.target = targetFH 
     """Output directory should be _nlin """
     self.nlinDir = nlinOutputDir
     """Empty array that we will fill with averages as we create them"""
     self.nlinAverages = [] 
     """Create the blurring resolution from the file resolution"""
     try: # the attempt to access the minc volume will fail if it doesn't yet exist at pipeline creation
         self.fileRes = rf.getFinestResolution(self.target)
     except: 
         # if it indeed failed, get resolution from the original file specified for 
         # one of the input files, which should exist. 
         # Can be overwritten by the user through specifying a nonlinear protocol.
         self.fileRes = rf.getFinestResolution(self.inputs[0].inputFileName)
     
     """
         Set default parameters before checking to see if a non-linear protocol has been
         specified. This is done first, since a non-linear protocol may specify only 
         a subset of the parameters, but all parameters must be set for the registration
         to run properly. 
         
         After default parameters are set, check for a specified non-linear protocol and
         override these parameters accordingly. Currently, this protocol must
         be a csv file that uses a SEMI-COLON to separate the fields. Examples are:
         pydpiper_apps_testing/test_data/minctracc_example_protocol.csv
         pydpiper_apps_testing/test_data/mincANTS_example_protocol.csv
        
         Each row in the csv is a different input to the either a minctracc or mincANTS call
         Although the number of entries in each row (e.g. generations) is variable, the
         specific parameters are fixed. For example, one could specify a subset of the
         allowed parameters (e.g. blurs only) but could not rename any parameters
         or use additional ones that haven't already been defined without subclassing. See
         documentation for additional details. 
         
         Note that if no protocol is specified, then defaults will be used. 
         Based on the length of these parameter arrays, the number of generations is set. 
     """
     
     self.defaultParams()
     if nlin_protocol:
         self.setParams(nlin_protocol)
     self.generations = self.getGenerations()   
     
     # Create new nlin group for each input prior to registration
     for i in range(len(self.inputs)):
         self.inputs[i].newGroup(groupName="nlin")

def MAGeTMask(atlases, inputs, numAtlases, regMethod, lsq12_protocol=None, nlin_protocol=None):
    """ Masking algorithm is as follows:
        1. Run HierarchicalMinctracc or mincANTS with mask=True, 
           using masks instead of labels. 
        2. Do voxel voting to find the best mask. (Or, if single atlas,
            use that transform)
        3. mincMath to multiply original input by mask to get _masked.mnc file
            (This is done for both atlases and inputs, though for atlases, voxel
             voting is not required.)
        4. Replace lastBasevol with masked version, since once we have created
            mask, we no longer care about unmasked version. 
        5. Clear out labels arrays, which were used to keep track of masks,
            as we want to re-set them for actual labels.
                
        Note: All data will be placed in a newly created masking directory
        to keep it separate from data generated during actual MAGeT. 
        """
    p = Pipeline()
    for atlasFH in atlases:
        maskDirectoryStructure(atlasFH, masking=True)
    for inputFH in inputs:
        maskDirectoryStructure(inputFH, masking=True)
        for atlasFH in atlases:
            sp = MAGeTRegister(inputFH, 
                               atlasFH, 
                               regMethod, 
                               name="initial", 
                               createMask=True,
                               lsq12_protocol=lsq12_protocol,
                               nlin_protocol=nlin_protocol)
            p.addPipeline(sp)          
    """ Prior to final masking, set log and tmp directories as they were."""
    for atlasFH in atlases:
        """Retrieve labels for use in new group. Assume only one"""
        labels = atlasFH.returnLabels(True)
        maskDirectoryStructure(atlasFH, masking=False)
        mp = maskFiles(atlasFH, True)
        p.addPipeline(mp)
        atlasFH.newGroup()
        atlasFH.addLabels(labels[0], inputLabel=True)
    for inputFH in inputs:
        maskDirectoryStructure(inputFH, masking=False)
        mp = maskFiles(inputFH, False, numAtlases)
        p.addPipeline(mp)
        # this will remove the "inputLabels"; labels that
        # come directly from the atlas library
        inputFH.clearLabels(True)
        # this will remove the "labels"; second generation
        # labels. I.e. labels from labels from the atlas library
        inputFH.clearLabels(False) 
        inputFH.newGroup()  
    return(p)    

def maskFiles(FH, isAtlas, numAtlases=1):
    """ Assume that if there is more than one atlas, multiple
        masks were generated and we need to perform a voxel_vote. 
        Otherwise, assume we are using inputLabels from crossing with
        only one atlas. 
    """
    #MF TODO: Make this more general to handle pairwise option. 
    p = Pipeline()
    if not isAtlas:
        if numAtlases > 1:
            voxel = voxelVote(FH, False, True)
            p.addStage(voxel)
            mincMathInput = voxel.outputFiles[0]  
        else:
            mincMathInput = FH.returnLabels(True)[0]
        FH.setMask(mincMathInput)
    else:
        mincMathInput = FH.getMask()
    mincMathOutput = fh.createBaseName(FH.resampledDir, FH.basename)
    mincMathOutput += "_masked.mnc"   
    logFile = fh.logFromFile(FH.logDir, mincMathOutput)
    cmd = ["mincmath"] + ["-clobber"] + ["-mult"]
    # In response to issue #135
    # the order of the input files to mincmath matters. By default the
    # first input files is used as a "like file" for the output file. 
    # We should make sure that the mask is not used for that, because
    # it has an image range from 0 to 1; not something we want to be
    # set for the masked output file
    #            average                              mask
    cmd += [InputFile(FH.getLastBasevol())] + [InputFile(mincMathInput)]
    cmd += [OutputFile(mincMathOutput)]
    mincMath = CmdStage(cmd)
    mincMath.setLogFile(LogFile(logFile))
    p.addStage(mincMath)
    FH.setLastBasevol(mincMathOutput)
    return(p)

 def __init__(self, inputFH, targetFH, statsKernels, additionalXfm=None):
     self.p = Pipeline()
     self.inputFH = inputFH
     self.targetFH = targetFH
     self.blurs = []
     self.setupBlurs(statsKernels)
     self.statsGroup = StatsGroup()
     self.setupXfms()
     """ additionalXfm is an optional transform that may be specified. If it is, 
         it is concatenated with the lastXfm from input to target. This additional
         transform must also be in the same direction as the lastXfm (e.g. input to target)
         Example usage: if the lastXfm from input to target goes from lsq12 to nlin space
         and you would like to calculate the absolute jacobians to lsq6 space, the additional
         transform specified may be the lsq6 to lsq12 transform from input to target. 
     """
     self.additionalXfm = additionalXfm
     self.fullStatsCalc()

 def __init__(self, inputFH, targetFH, blurs, inputArray=None, scalingFactor=None):
     self.p = Pipeline()
     self.inputFH = inputFH
     self.targetFH = targetFH
     self.blurs = blurs
     self.statsGroup = StatsGroup()
     self.setupXfms()
     """Optional inputArray used to calculate an average displacement and use for recentering."""
     if inputArray:
         self.dispToAvg = []
         self.setupDispArray(inputArray)
     else:
         self.dispToAvg = None
     """   
         Specify an optional xfm to be used when calculating the 
         scaled jacobians. This jacobian will then be concatenated with the
         self.linearXfm, the linear portion of the final non-linear transform from input to target.
         
         A  
         
     """
     self.scalingFactor = scalingFactor

class AbstractApplication(object):
    """Framework class for writing applications for PydPiper. 
    
       This class defines the default behaviour for accepting common command-line options, and executing the application
       under various queueing systems. 
       
       Subclasses should extend the following methods:
           setup_appName()
           setup_logger() [optional, default method is defined here]
           setup_options()
           run()
    
       Usage: 
          class MyApplication(AbstractApplication):
                ... 
           
          if __name__ == "__main__":
              application = MyApplication()
              application.start()
    """
    def __init__(self):
        # use an environment variable to look for a default config file
        # Alternately, we could use a default location for the file
        # (say `files = ['/etc/pydpiper.cfg', '~/pydpiper.cfg', './pydpiper.cfg']`)
        default_config_file = os.getenv("PYDPIPER_CONFIG_FILE")
        if default_config_file is not None:
            files = [default_config_file]
        else:
            files = []
        self.parser = MyParser(default_config_files=files)
        self.__version__ = get_distribution("pydpiper").version  # pylint: disable=E1101
    
    def _setup_options(self):
            # PydPiper options
        addExecutorArgumentGroup(self.parser)
        addApplicationArgumentGroup(self.parser)
    
    def _print_version(self):
        if self.options.show_version:
            print(self.__version__)
            sys.exit()
    
    def _setup_pipeline(self, options):
        self.pipeline = Pipeline(options)

    # FIXME check that only one server is running with a given output directory
    def _setup_directories(self):
        """Output and backup directories setup here."""
        if not self.options.output_directory:
            self.outputDir = os.getcwd()
        else:
            self.outputDir = makedirsIgnoreExisting(self.options.output_directory)
        self.pipeline.setBackupFileLocation(self.outputDir)

    def reconstructCommand(self):    
        reconstruct = ' '.join(sys.argv)
        logger.info("Command is: " + reconstruct)
        logger.info("Command version : " + self.__version__)
        # also, because this is probably a better file for it (also has similar
        # naming conventions as the pipeline-stages.txt file:
        fileForCommandAndVersion = self.options.pipeline_name + "-command-and-version-" + time.strftime("%d-%m-%Y-at-%H-%m-%S") + ".sh"
        pf = open(fileForCommandAndVersion, "w")
        pf.write("#!/usr/bin/env bash\n")
        pf.write("# Command version is: " + self.__version__ + "\n")
        pf.write("# Command was: \n")
        pf.write(reconstruct + '\n')
        pf.write("# options were: \n# %s" % self.options)
        pf.close()
        
    def start(self):
        logger.info("Calling `start`")
        self._setup_options()
        self.setup_options()
        self.options = self.parser.parse_args()
        self.args = self.options.files

        self._print_version()

        # Check to make sure some executors have been specified if we are 
        # actually going to run:
        if self.options.execute:
            noExecSpecified(self.options.num_exec)
             
        self._setup_pipeline(self.options)
        self._setup_directories()
        
        self.appName = self.setup_appName()
        self.setup_logger()

        # TODO this doesn't capture environment variables
        # or contents of any config file so isn't really complete
        self.reconstructCommand()

        pbs_submit = self.options.queue_type == "pbs" \
                     and not self.options.local

        # --create-graph causes the pipeline to be constructed
        # both at PBS submit time and on the grid; this may be an extremely
        # expensive duplication
        if (self.options.execute and not pbs_submit) or self.options.create_graph:
#.........这里部分代码省略.........

 def _setup_pipeline(self, options):
     self.pipeline = Pipeline(options)

class CalcStats(object):
    """Statistics calculation between an input and target. 
       This class calculates multiple displacement fields, relative and absolute jacobians.  
       General functionality as follows:
       1. Class instantiated with input, target and statsKernels. Note that here, the statsKernels 
          specified are blurs used to smooth the displacement fields prior to additional calculations. 
          They may be a string of comma separated values or an array of doubles. 
       2. An additional transform may also be included to calculate absolute 
          jacobians to a different space, as is described in the __init__ function, 
          documentation and elsewhere in the code.  
       3. If needed, invert transform between input and target in setupXfms(). This is necessary
          as this class assumes that the target is the reference space, from which all stats
          are calculated. 
       4. Call fullStatsCalc. This calculates linear and 
          pure nonlinear displacement before calculating jacobians.    
       5. Ability to recenter displacements using an average may be re-added in the future. 
    """
    def __init__(self, inputFH, targetFH, statsKernels, additionalXfm=None):
        self.p = Pipeline()
        self.inputFH = inputFH
        self.targetFH = targetFH
        self.blurs = []
        self.setupBlurs(statsKernels)
        self.statsGroup = StatsGroup()
        self.setupXfms()
        """ additionalXfm is an optional transform that may be specified. If it is, 
            it is concatenated with the lastXfm from input to target. This additional
            transform must also be in the same direction as the lastXfm (e.g. input to target)
            Example usage: if the lastXfm from input to target goes from lsq12 to nlin space
            and you would like to calculate the absolute jacobians to lsq6 space, the additional
            transform specified may be the lsq6 to lsq12 transform from input to target. 
        """
        self.additionalXfm = additionalXfm
        self.fullStatsCalc()
        
    def setupBlurs(self, statsKernels):
        if isinstance(statsKernels, list):
            self.blurs = statsKernels
        elif isinstance(statsKernels, str):
            for i in statsKernels.split(","):
                self.blurs.append(float(i))
        else:
            print("Improper type of blurring kernels specified for stats calculation: " + str(statsKernels))
            sys.exit()
           
    def setupXfms(self):
        self.xfm = self.inputFH.getLastXfm(self.targetFH)
        if not self.xfm:
            print("Cannot calculate statistics. No transform between input and target specified.")
            print("Input: " + self.inputFH.getLastBasevol())
            print("Target: " + self.targetFH.getLastBasevol())
            sys.exit()
        else:
            self.invXfm = self.targetFH.getLastXfm(self.inputFH)
            if not self.invXfm:
                xi = xfmInvert(self.xfm, FH=self.inputFH)
                self.p.addStage(xi)
                self.invXfm = xi.outputFiles[0]
    
    def fullStatsCalc(self):
        self.linAndNlinDisplacement()
        self.calcDetAndLogDet(useFullDisp=False)  # Calculate relative jacobians
        self.calcDetAndLogDet(useFullDisp=True)   # Calculate absolute jacobians
    
    def calcFullDisplacement(self):
        """Calculate full displacement from target to input. If an
           additionaXfm is specified, it is concatenated to self.xfm here """
        if self.additionalXfm:
            outXfm = createOutputFileName(self.inputFH, self.xfm, "transforms", "_with_additional.xfm")
            xc = xfmConcat([self.additionalXfm, self.xfm], outXfm, fh.logFromFile(self.inputFH.logDir, outXfm))
            self.p.addStage(xc)
            xi = xfmInvert(xc.outputFiles[0], FH=self.inputFH)
            self.p.addStage(xi)
            fullDisp = mincDisplacement(self.targetFH, self.inputFH, transform=xi.outputFiles[0])
        else:
            fullDisp = mincDisplacement(self.targetFH, self.inputFH, transform=self.invXfm)
        self.p.addStage(fullDisp)
        self.fullDisp = fullDisp.outputFiles[0]
    
    def calcNlinDisplacement(self):
        """Calculate pure non-linear displacement from target to input
           1. Concatenate self.invXfm (target to input xfm) and self.linearPartOfNlinXfm   
           2. Compute mincDisplacement on this transform. 
        """
        pureNlinXfm = createOutputFileName(self.inputFH, self.invXfm, &qu