def nlin_part(xfm : XfmHandler, inv_xfm : Optional[XfmHandler] = None) -> Result[XfmHandler]:
    """
    *** = non linear deformations
    --- = linear (affine) deformations

    Input:
    xfm     :     ******------>
    inv_xfm :    <******------ (optional)

    Calculated:
    inv_lin_xfm :      <------

    Returned:
    concat :      ******------> +
                       <------
    equals :      ******>

    Compute the nonlinear part of a transform as follows:
    go forwards across xfm and then backwards across the linear part
    of the inverse xfm (by first calculating the inverse or using the one supplied) 
    Finally, use minc_displacement to compute the resulting gridfile of the purely 
    nonlinear part.

    The optional inv_xfm (which must be the inverse!) is an optimization -
    we don't go looking for an inverse by filename munging and don't programmatically
    keep a log of operations applied, so any preexisting inverse must be supplied explicitly.
    """
    s = Stages()
    inv_xfm = inv_xfm or s.defer(invert_xfmhandler(xfm))
    inv_lin_part = s.defer(lin_from_nlin(inv_xfm)) 
    xfm = s.defer(concat_xfmhandlers([xfm, inv_lin_part]))
    return Result(stages=s, output=xfm)

def mincblob(op : str, grid : MincAtom, subdir : str = "tmp") -> Result[MincAtom]:
    """
    Low-level mincblob wrapper with the one exception being the determinant option. By
    default the inner clockwork of mincblob subtracts 1 from all determinant values that
    are being calculated. As such, 1 needs to be added to the result of the mincblob call.
    We will do that here, because it makes most sense here.
    >>> stages = mincblob('determinant', MincAtom("/images/img_1.mnc", pipeline_sub_dir="/tmp")).stages
    >>> [s.render() for s in stages]
    ['mincblob -clobber -determinant /images/img_1.mnc /tmp/img_1/img_1_determinant.mnc']
    """
    if op not in ["determinant", "trace", "translation", "magnitude"]:
        raise ValueError('mincblob: invalid operation %s' % op)

    # if we are calculating the determinant, the first file produced is a temp file:
    if op == "determinant":
        out_file = grid.newname_with_suffix("_temp_det", subdir=subdir)
    else:
        out_file = grid.newname_with_suffix('_' + op, subdir=subdir)

    stage = CmdStage(inputs=(grid,), outputs=(out_file,),
                 cmd=['mincblob', '-clobber', '-' + op, grid.path, out_file.path])

    s = Stages([stage])
    # now create the proper determinant if that's what was asked for
    if op == "determinant":
        result_file = s.defer(mincmath(op='add',
                                       const=1,
                                       vols=[out_file],
                                       subdir=subdir,
                                       new_name=grid.filename_wo_ext + "_det"))
    else:
        result_file = out_file

    return Result(stages=s, output=result_file)

def lsq6_pipeline(options):
    # TODO could also allow pluggable pipeline parts e.g. LSQ6 could be substituted out for the modified LSQ6
    # for the kidney tips, etc...
    output_dir    = options.application.output_directory
    pipeline_name = options.application.pipeline_name

    # TODO this is tedious and annoyingly similar to the registration chain and MBM ...
    lsq6_dir      = os.path.join(output_dir, pipeline_name + "_lsq6")
    processed_dir = os.path.join(output_dir, pipeline_name + "_processed")
    imgs = get_imgs(options.application)

    s = Stages()

    # TODO this is quite tedious and duplicates stuff in the registration chain ...
    resolution = (options.registration.resolution or
                  get_resolution_from_file(
                      s.defer(registration_targets(lsq6_conf=options.lsq6,
                                                   app_conf=options.application,
                                                   reg_conf=options.registration)).registration_standard.path))

    # FIXME: why do we have to call registration_targets *outside* of lsq6_nuc_inorm? is it just because of the extra
    # options required?
    targets = s.defer(registration_targets(lsq6_conf=options.lsq6,
                                   app_conf=options.application,
                                   reg_conf=options.registration,
                                   first_input_file=imgs[0]))
    # This must happen after calling registration_targets otherwise it will resample to options.registration.resolution
    options.registration = options.registration.replace(resolution=resolution)
    lsq6_result = s.defer(lsq6_nuc_inorm(imgs=imgs,
                                         resolution=resolution,
                                         registration_targets=targets,
                                         lsq6_dir=lsq6_dir,
                                         lsq6_options=options.lsq6))

    return Result(stages=s, output=lsq6_result)

def asymmetry_pipeline(options):

    output_dir    = options.application.output_directory
    pipeline_name = options.application.pipeline_name
    processed_dir = os.path.join(output_dir, pipeline_name + "_processed")

    s = Stages()

    #imgs_ = [MincAtom(f, pipeline_sub_dir=processed_dir) for f in options.application.files]

    imgs_ = get_imgs(options.application)

    check_MINC_input_files([img.path for img in imgs_])

    imgs  = pd.Series(imgs_, index=[img.filename_wo_ext for img in imgs_])
    flipped_imgs = imgs.apply(lambda img: s.defer(volflip(img)))  # TODO add flags to control flip axis ...

    # TODO ugly - MincAtom API should allow this somehow without mutation (also, how to pass into `volflip`, etc.?)
    for f_i in flipped_imgs:
        f_i.output_sub_dir += "_flipped"

    check_MINC_input_files(imgs.apply(lambda img: img.path))

    grouped_files_df = pd.DataFrame({'file' : pd.concat([imgs, flipped_imgs])}).assign(group=lambda df: df.index)

    two_level_result = s.defer(two_level(grouped_files_df, options=options))

    return Result(stages=s, output=two_level_result)

def cortical_thickness_pipeline(options):
    s = Stages()

    #imgs = [MincAtom(name, pipeline_sub_dir=os.path.join(options.application.output_directory,
    #                                                     options.application.pipeline_name + "_processed"))
    #        for name in options.application.files]

    pipeline_sub_dir = os.path.join(options.application.output_directory,
                                    options.application.pipeline_name + "_processed")

    #def atom(atom_type, file):
    #    return atom_type(file, pipeline_sub_dir=pipeline_sub_dir)  # TODO output_sub_dir, ....

    # TODO are all these fields actually used?  If not, omit from CSV?
    xfms = (pd.read_csv(options.thickness.xfm_csv)
            .apply(axis=1,  # TODO fill out <..>Atom(...) fields ...
                   func=lambda row: XfmHandler(
                          source=MincAtom(row.source, pipeline_sub_dir=pipeline_sub_dir),
                          target=MincAtom(row.target, pipeline_sub_dir=pipeline_sub_dir),
                          resampled=None,   #MincAtom(row.resampled, pipeline_sub_dir=pipeline_sub_dir),
                          xfm=XfmAtom(row.xfm, pipeline_sub_dir=pipeline_sub_dir))))
    # TODO better way to unpack?

    result = s.defer(cortical_thickness(xfms=xfms,
                                        atlas=NotImplemented,
                                        label_mapping=options.thickness.label_mapping,
                                        atlas_fwhm=options.thickness.atlas_fwhm,
                                        thickness_fwhm=options.thickness.thickness_fwhm))

    return Result(stages=s, output=result)

def tamarack_pipeline(options):

    output_dir    = options.application.output_directory
    pipeline_name = options.application.pipeline_name
    #processed_dir = os.path.join(output_dir, pipeline_name + "_processed")
    first_level_dir = os.path.join(output_dir, pipeline_name + "_first_level")

    s = Stages()

    with open(options.application.csv_file, 'r') as f:
        files_df = (pd.read_csv(filepath_or_buffer=f,
                                usecols=['group', 'filename'])
                    .assign(file=lambda df:
                                   df.apply(axis="columns",
                                            func=lambda r:
                                                   MincAtom(r.filename.strip(),
                                                            pipeline_sub_dir=os.path.join(first_level_dir,
                                                                                          "%s_processed" % r.group.strip())))))

    check_MINC_input_files(files_df.file.apply(lambda img: img.path))

    #grouped_files_df = pd.DataFrame({'file' : pd.concat([imgs])}).assign(group=lambda df: df.index)

    tamarack_result = s.defer(tamarack(files_df, options=options))

    tamarack_result.first_level_results.applymap(maybe_deref_path).to_csv("first_level_results.csv", index=False)
    tamarack_result.resampled_determinants.applymap(maybe_deref_path).to_csv("resampled_determinants.csv", index=False)
    tamarack_result.overall_determinants.applymap(maybe_deref_path).to_csv("overall_determinants.csv", index=False)

    return Result(stages=s, output=tamarack_result)

 def scale_transform(xfm, scale, newname_wo_ext):
     s = Stages()
     defs = s.defer(as_deformation(transform=xfm.xfm, reference=xfm.source))
     scaled_defs = (defs.xfm.newname(newname_wo_ext) if newname_wo_ext else
                     defs.xfm.newname_with_suffix("_scaled_%s" % scale))
     s.defer(CmdStage(cmd=['c3d', '-scale', str(scale), defs.path, "-o", scaled_defs.path],
                      inputs=(defs,), outputs=(scaled_defs,)))
     return Result(stages=s, output=scaled_defs)

def determinant(displacement_grid : MincAtom) -> Result[MincAtom]:
    """
    Takes a displacement field (deformation grid, vector field, those are
    all the same thing) and calculates the proper determinant (mincblob()
    takes care of adding 1 to the silly output of running mincblob directly)
    """
    s = Stages()
    det = s.defer(mincblob(op='determinant', grid=displacement_grid))
    return Result(stages=s, output=det)

def convert(infile : ImgAtom, out_ext : str) -> Result[ImgAtom]:
    s = Stages()
    outfile = infile.newext(ext=out_ext)
    if infile.mask is not None:
        outfile.mask = s.defer(convert(infile.mask, out_ext=out_ext))
    if infile.labels is not None:
        outfile.mask = s.defer(convert(infile.labels, out_ext=out_ext))
    s.add(CmdStage(inputs=(infile,), outputs=(outfile,),
                   cmd = ['c3d', infile.path, '-o', outfile.path]))
    return Result(stages=s, output=outfile)

 def to_mni_xfm(xfm):
     s = Stages()
     defs = xfm.newname_with_suffix("_defs", subdir="tmp")
     s.add(CmdStage(cmd=["transformix", "-def", "all",
                         "-out", defs.dir,
                         "-tp", xfm.path,
                         "-xfm", os.path.join(defs.filename_wo_ext, defs.ext)],
                    inputs=(xfm,), outputs=(defs,)))
     out_xfm = s.defer(itk.itk_convert_xfm(defs, out_ext=".mnc"))
     return Result(stages=s, output=out_xfm)

 def average_transforms(xfms, avg_xfm):
     s = Stages()
     defs = [s.defer(as_deformation(transform=xfm.xfm, reference_image=xfm.source)) for xfm in xfms]
     #avg_img = NotImplemented
     avg = imageToXfm(s.defer(average_images(defs,
                                             avg_file=xfmToImage(avg_xfm),
                                             #output_dir=os.path.join(defs[0].pipeline_sub_dir,
                                             #                        defs[0].output_sub_dir,
                                             #                        "transforms")
                                             )))
     return Result(stages=s, output=avg)

def NLIN_pipeline(options):

    # if options.application.files is None:
    #     raise ValueError("Please, some files! (or try '--help')")  # TODO make a util procedure for this

    output_dir    = options.application.output_directory
    pipeline_name = options.application.pipeline_name

    # TODO this is tedious and annoyingly similar to the registration chain and MBM and LSQ6 ...
    processed_dir = os.path.join(output_dir, pipeline_name + "_processed")
    nlin_dir      = os.path.join(output_dir, pipeline_name + "_nlin")

    resolution = (options.registration.resolution  # TODO does using the finest resolution here make sense?
                  or min([get_resolution_from_file(f) for f in options.application.files]))

    imgs = get_imgs(options.application)

    # imgs = [MincAtom(f, pipeline_sub_dir=processed_dir) for f in options.application.files]

    # determine NLIN settings by overriding defaults with
    # any settings present in protocol file, if it exists
    # could add a hook to print a message announcing completion, output files,
    # add more stages here to make a CSV

    initial_target_mask = MincAtom(options.nlin.target_mask) if options.nlin.target_mask else None
    initial_target = MincAtom(options.nlin.target, mask=initial_target_mask)

    full_hierarchy = get_nonlinear_configuration_from_options(nlin_protocol=options.nlin.nlin_protocol,
                                                              reg_method=options.nlin.reg_method,
                                                              file_resolution=resolution)

    s = Stages()

    nlin_result = s.defer(nlin_build_model(imgs, initial_target=initial_target, conf=full_hierarchy, nlin_dir=nlin_dir))

    # TODO return these?
    inverted_xfms = [s.defer(invert_xfmhandler(xfm)) for xfm in nlin_result.output]

    if options.stats.calc_stats:
        # TODO: put the stats part behind a flag ...

        determinants = [s.defer(determinants_at_fwhms(
                                  xfm=inv_xfm,
                                  inv_xfm=xfm,
                                  blur_fwhms=options.stats.stats_kernels))
                        for xfm, inv_xfm in zip(nlin_result.output, inverted_xfms)]

        return Result(stages=s,
                      output=Namespace(nlin_xfms=nlin_result,
                                       avg_img=nlin_result.avg_img,
                                       determinants=determinants))
    else:
        # there's no consistency in what gets returned, yikes ...
        return Result(stages=s, output=Namespace(nlin_xfms=nlin_result, avg_img=nlin_result.avg_img))

 def f(imgs: List[MincAtom],
       initial_target: MincAtom,
       conf: reg_module.Conf,
       nlin_dir: str,
       nlin_prefix: str,
       tournament_name_wo_ext: str = "tournament") -> Result[List[XfmHandler]]:
   s = Stages()
   Weight = int
   def h(xfms : List[XfmHandler], name_wo_ext : str) -> List[XfmHandler]:
       # TODO add weights to each return
       # TODO check len(...) == 0 case??
       if len(xfms) <= 1:
           return xfms
       else:
           first_half  = xfms[: len(xfms)//2]
           second_half = xfms[len(xfms)//2 :]
           first_half_result  = h(first_half,  name_wo_ext=name_wo_ext + "_L")
           second_half_result = h(second_half, name_wo_ext=name_wo_ext + "_R")
           A_halfway_to_B, B_halfway_to_A, avg_img = s.defer(
               nonlinear_midpoint_xfm(
                   img_A = first_half_result[0].resampled,
                   img_B = second_half_result[0].resampled,
                   out_name_wo_ext=name_wo_ext,
                   nlin_algorithm=reg_module,
                   conf=conf,
                   out_dir=nlin_dir))
           xfms_to_midpoint = ([XfmHandler(source=xfm.source,
                                           target=avg_img,
                                           resampled=A_halfway_to_B.resampled,
                                           xfm=s.defer(xfmconcat([xfm.xfm, A_halfway_to_B.xfm],
                                                                 name="%s_%s" % (xfm.source.filename_wo_ext,
                                                                                 name_wo_ext))))
                                for xfm in first_half_result]
                               + [XfmHandler(source=xfm.source,
                                             target=avg_img,
                                             resampled=B_halfway_to_A.resampled,
                                             xfm=s.defer(xfmconcat([xfm.xfm, B_halfway_to_A.xfm],
                                                                   name="%s_%s" % (xfm.source.filename_wo_ext,
                                                                                   name_wo_ext))))
                                  for xfm in second_half_result])
           return xfms_to_midpoint
   identity_xfm = s.defer(param2xfm(out_xfm=XfmAtom(pipeline_sub_dir=imgs[0].pipeline_sub_dir,
                                                    output_sub_dir=imgs[0].output_sub_dir,
                                                    name=os.path.join(imgs[0].pipeline_sub_dir,
                                                                      imgs[0].output_sub_dir,
                                                                      "id.xfm"))))
   initial_xfms = [XfmHandler(source=img, target=img,
                              resampled=img, xfm=identity_xfm) for img in imgs]
   xfms_to_avg = h(initial_xfms, tournament_name_wo_ext)
   avg_img = xfms_to_avg[0].target
   return Result(stages=s, output=WithAvgImgs(avg_img=avg_img, avg_imgs=[avg_img],
                                              output=xfms_to_avg))

 def average_transforms(xfms, avg_xfm):
     intermediate_xfm = avg_xfm.newname_with_suffix("_inter", subdir="tmp")
     s = Stages()
     s.add(CmdStage(cmd=["echo", ('(Transform "WeightedCombinationTransform")\n'
                                  '(SubTransforms %s)\n'
                                  '(NormalizeCombinationsWeights "true")\n') %
                                    ' '.join(sorted(xfm.path for xfm in xfms))],
                    inputs=xfms, outputs=(intermediate_xfm,)))
     s.add(CmdStage(cmd=["transformix", "-def", "all",
                         "-out", os.path.dirname(avg_xfm.path),
                         "-tp", intermediate_xfm.path,
                         "-xfm", avg_xfm.path],
                    inputs=(intermediate_xfm,), outputs=(avg_xfm,)))

def nlin_displacement(xfm : XfmHandler, inv_xfm : Optional[XfmHandler] = None) -> Result[MincAtom]:
    """
    See: nlin_part().

    This returns the nonlinear part of the input
    transformation (xfm) in the form of a grid file (vector field).
    All transformations are encapsulated in this field (linear parts
    that are normally specified in the .xfm file are placed in the
    vector field)
    """
    
    s = Stages()
    return Result(stages=s,
                  output=s.defer(minc_displacement(
                                   s.defer(nlin_part(xfm, inv_xfm=inv_xfm)))))

def resample(img,
             xfm,  # TODO: update to handler?
             like,
             invert = False,
             use_nn_interpolation = None,
             new_name_wo_ext: str = None,
             subdir: str = None,
             postfix: str = None):

    s = Stages()

    if not subdir:
        subdir = 'resampled'

    # we need to get the filename without extension here in case we have
    # masks/labels associated with the input file. When that's the case,
    # we supply its name with "_mask" and "_labels" for which we need
    # to know what the main file will be resampled as
    if not new_name_wo_ext:
        # FIXME this is wrong when invert=True
        new_name_wo_ext = xfm.filename_wo_ext + '-resampled'

    new_img = s.defer(resample_simple(img=img, xfm=xfm, like=like,
                                      invert=invert,
                                      use_nn_interpolation=use_nn_interpolation,
                                      new_name_wo_ext=new_name_wo_ext,
                                      subdir=subdir))
    new_img.mask = s.defer(resample_simple(img=img.mask, xfm=xfm, like=like,
                                           use_nn_interpolation=True,
                                           invert=invert,
                                           new_name_wo_ext=new_name_wo_ext + "_mask",
                                           subdir=subdir)) if img.mask is not None else None
    new_img.labels = s.defer(resample_simple(img=img.labels, xfm=xfm, like=like,
                                             use_nn_interpolation=True,
                                             invert=invert,
                                             new_name_wo_ext=new_name_wo_ext + "_labels",
                                             subdir=subdir)) if img.labels is not None else None

    # Note that new_img can't be used for anything until the mask/label files are also resampled.
    # This shouldn't create a problem with stage dependencies as long as masks/labels appear in inputs/outputs of CmdStages.
    # (If this isn't automatic, a relevant helper function would be trivial.)
    # TODO: can/should this be done semi-automatically? probably ...
    return Result(stages=s, output=new_img)

def average_images(imgs        : Sequence[ImgAtom],
                   dimensions  : int = 3,
                   normalize   : bool = False,
                   output_dir  : str = '.',
                   name_wo_ext : str = "average",
                   out_ext     : Optional[str] = None,
                   avg_file    : Optional[ITKImgAtom] = None) -> Result[ITKImgAtom]:

    s = Stages()

    if len(imgs) == 0:
        raise ValueError("`AverageImages` arg `imgs` is empty (can't average zero files)")

    ext = out_ext or imgs[0].ext

    # the output_dir basically gives us the equivalent of the pipeline_sub_dir for
    # regular input files to a pipeline, so use that here
    avg = avg_file or ImgAtom(name=os.path.join(output_dir, '%s.todo' % name_wo_ext),
                              orig_name=None,
                              pipeline_sub_dir=output_dir)
    avg.ext = ext

    # if all input files have masks associated with them, add the combined mask to
    # the average:
    # TODO what if avg_file has a mask ... should that be used instead? (then rename avg -> avg_file above)
    all_inputs_have_masks = all((img.mask for img in imgs))
    if all_inputs_have_masks:
        combined_mask = (ImgAtom(name=os.path.join(avg_file.dir, '%s_mask.todo' % avg_file.filename_wo_ext),
                                 orig_name=None,
                                 pipeline_sub_dir=avg_file.pipeline_sub_dir)
                         if avg_file is not None else
                         ImgAtom(name=os.path.join(output_dir, '%s_mask.todo' % name_wo_ext),
                                 orig_name=None,
                                 pipeline_sub_dir=output_dir))
        combined_mask.ext = ext
        s.defer(max(imgs=sorted({img_inst.mask for img_inst in imgs}),
                    out_img=combined_mask))
        avg.mask = combined_mask
    s.add(CmdStage(inputs = imgs,
                   outputs = (avg,),
                   cmd = ["AverageImages", str(dimensions), avg.path, "%d" % normalize]
                         + [img.path for img in imgs]))
    return Result(stages=s, output=avg)

def stage_embryos_pipeline(options):
    s = Stages()

    imgs = get_imgs(options.application)
    rough_volume_imgs = get_volume_estimate(imgs)
    imgs_and_rough_volume = pd.DataFrame({"mincatom" : imgs,
                                          "rough_volume" : pd.Series(rough_volume_imgs, dtype=float)})

    check_MINC_input_files([img.path for img in imgs])

    output_directory = options.application.output_directory
    output_sub_dir = os.path.join(output_directory,
                                  options.application.pipeline_name + "_4D_atlas")

    time_points_in_4D_atlas = instances_in_4D_atlas_from_csv(options.staging.staging.csv_4D,
                                                             output_sub_dir)

    # we can use the resolution of one of the time points in the 4D atlas
    # for all the registrations that will be run.
    resolution = get_resolution_from_file(time_points_in_4D_atlas.loc[0]["mincatom"].orig_path)

    print(options.staging.lsq12)

    lsq12_conf = get_linear_configuration_from_options(options.staging.lsq12,
                                                       transform_type=LinearTransType.lsq12,
                                                       file_resolution=resolution)

    nlin_component = get_nonlinear_component(options.staging.nlin.reg_method)

    # match each of the embryos individually
    for i in range(imgs_and_rough_volume.shape[0]):
        s.defer(match_embryo_to_4D_atlas(imgs_and_rough_volume.loc[i],
                                         time_points_in_4D_atlas,
                                         lsq6_conf=options.staging.lsq6,
                                         lsq12_conf=lsq12_conf,
                                         nlin_module=nlin_component,
                                         resolution=resolution,
                                         nlin_options=options.staging.nlin))


    return Result(stages=s, output=None)

        def build_model(imgs,
                        conf,
                        nlin_dir,
                        nlin_prefix,
                        initial_target,
                        output_name_wo_ext = None):
            s = Stages()
            mincify = base_build_model.ToMinc
            imgs = tuple(s.defer(mincify.from_mnc(img)) for img in imgs)
            result = s.defer(base_build_model.build_model(imgs=imgs, conf=conf,
                                                  nlin_dir=nlin_dir, nlin_prefix=nlin_prefix,
                                                  initial_target=s.defer(mincify.from_mnc(initial_target))
                                                  #output_name_wo_ext=output_name_wo_ext
                                                  ))

            def wrap_output_xfmh(xfmh):
                return XfmHandler(source=s.defer(mincify.to_mnc(xfmh.source)) if xfmh.source else None,
                                  target=s.defer(mincify.to_mnc(xfmh.target)) if xfmh.target else None,
                                  resampled=s.defer(mincify.to_mnc(xfmh.resampled)) if xfmh.has_resampled() else None,
                                  xfm=s.defer(mincify.to_mni_xfm(xfmh.xfm)),
                                  inverse=wrap_output_xfmh(xfmh.inverse) if xfmh.has_inverse() else None)

            return Result(stages=s, output=WithAvgImgs(avg_imgs=[s.defer(mincify.to_mnc(img))
                                                                 for img in result.avg_imgs],
                                                       avg_img=s.defer(mincify.to_mnc(result.avg_img)),
                                                       output=[wrap_output_xfmh(x)
                                                               for x in result.output]))

def match_embryo_to_4D_atlas(embryo_with_volume_est,
                             full_4D_atlas_info,
                             lsq6_conf: LSQ6Conf,
                             lsq12_conf: MinctraccConf,
                             nlin_module: NLIN,
                             resolution: float,
                             nlin_options):
    s = Stages()

    # 1 what's the closest match in the 4D atlas?
    mid_index = get_index_closest_volume_match(embryo_with_volume_est["rough_volume"].astype(float), full_4D_atlas_info)

    print("Best initial match for: \n", embryo_with_volume_est["mincatom"].orig_path, " ", full_4D_atlas_info.loc[mid_index]["timepoint"])

    # register embryo to closest match +/- 5 time points
    # make sure we don't index outside the possible range
    lowest_index  = max(0, mid_index - 7)
    highest_index = min(full_4D_atlas_info.shape[0] - 1, mid_index + 7)

    all_transforms = [s.defer(lsq6_lsq12_nlin(source=embryo_with_volume_est["mincatom"],
                                              target=full_4D_atlas_info.loc[i]["mincatom"],
                                              lsq6_conf=lsq6_conf,
                                              lsq12_conf=lsq12_conf,
                                              nlin_module=nlin_module,
                                              resolution=resolution,
                                              nlin_options=nlin_options.nlin_protocol,
                                              resampled_post_fix_string="E" + str(full_4D_atlas_info.loc[i]["timepoint"]))) for
                      i in range(lowest_index, highest_index + 1, 1)]

    # gather stats on those registrations
    # the match is determined by the sum of the magnitude
    # of the inverse transformation from 4D instance -> embryo
    # using the mask of the 4D instance to limit the total sum
    # 1) calculate inverse
    all_inv_transforms = [s.defer(invert_xfmhandler(xfm)) for xfm in all_transforms]
    minc_displacement_grids = [s.defer(minc_displacement(inv_xfm)) for inv_xfm in all_inv_transforms]
    magnitudes = [s.defer(mincblob(op='magnitude', grid=disp_grid)) for disp_grid in minc_displacement_grids]

    return Result(stages=s, output=all_transforms)