Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Parallelize conformation #666

Merged
merged 20 commits into from
Aug 21, 2017

Conversation

anibalsolon
Copy link
Contributor

Step t1_conform can become a bottleneck since image resampling is executed in series. However, it is possible to run it in parallel, for multiple T1w images.

This PR intends to parallelize this process by breaking it into two steps:

  1. Decide the zoom and scale that images should take, based on all images
  2. For each image, resample it to fit the new requirements

This process is executed before and after t1_merge, so a sub-workflow was created to wrap both steps.

This PR solves issue #613 .

@effigies
Copy link
Member

Hi, thanks for this. I think if we're going to do this at this point, we should probably simplify the t1_reorient step, rather than build a workflow. Most of the machinery is entirely skipped, because t1_reorient only ever gets one image as input. I'll have a look through this in a bit; I think we can manage this with pretty minimal changes to what you've added.

@anibalsolon
Copy link
Contributor Author

Hi @effigies ! Thank you for the feedback.

I believe it is possible to abstract the nb.as_closest_canonical as a Step, and put it before PruneExcessiveZoom, so the wf become:

  • MapNode "ClosestCanonical" (maybe better naming)
  • Node PruneExcessiveZoom
  • MapNode ConformSeries

and the t1_reorient do not use the whole wf, but just Node ClosestCanonical.

@effigies
Copy link
Member

I think the ClosestCanonical node would be perfectly reasonable, but there's no need to use it before and after the template creation; we can leave nb.as_closest_canonical in the PruneExcessiveZoom and ConformSeries interfaces; it's a basically free operation in memory, so there's no need to save the result. We can just run it after t1_template.

As to the name, maybe "Reorient", "OrientRAS"? I'm not too picky here.

Also, WRT "PruneExcessiveZoom", that's only one part of what that step does (and actually the rarest, in my experience). The main purpose it serves is to calculate the template dimensions. "TemplateDimensions", "CalcTemplateSize", "TemplateParameters", or something along those lines would be better.

What do you think?

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's a first pass review of the code. Once we decide how to go on these issues and see what the result looks like, I'll do a more fiddly, detailed review.

target_shape = traits.List(traits.Int,
desc='Target shape information')
target_span = traits.List(traits.Float,
desc='Target span information')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wouldn't bother with target_span, since it's just the dot-product of target_zooms and target_shape. Doing the multiplication avoids checking for the case where span != zooms * shape.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was in doubt to parametrize it or compute the dot product on the next step, so I will now compute it. 👌

t1w_valid_list = OutputMultiPath(exists=True, desc='valid T1w images')

target_zooms = traits.List(traits.Float,
desc='Target zoom information')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is a known length, you can either set a min/max length of 3, or do:

target_zooms = traits.Tuple(traits.Float, traits.Float, traits.Float, desc='...')

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perfect! I didn't know this format. Will use it now :)


class PruneExcessiveZoomOutputSpec(TraitedSpec):
t1w_dropped_list = OutputMultiPath(exists=True, desc='invalid T1w images')
t1w_valid_list = OutputMultiPath(exists=True, desc='valid T1w images')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The t1w_dropped_list isn't going to be used anywhere, is it? I'm not sure that this is going to be generally useful outside of fmriprep, so we don't need to add pieces we won't use. (If we move it into nipype because the people are clamoring for it, we can revisit.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes! Initially, it was to be used for a posterior reporting, but I ended adding reporting in this same interface. I'll remove it for now.


self._results['t1w_list'] = out_names
self._results['target_zooms'] = target_zooms.tolist()
self._results['target_shape'] = target_shape.astype(np.int64).tolist()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we need 64 bits to describe the shape of an array, we're in trouble. This should already be an int, so I think you can skip the astype.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hehe indeed! I cast it because it was using floats, but I will double check it, just in case. If it is not possible to remove this casting, I will switch np.int64 to int.

t1w = File(exists=True, desc='conformed T1w image')


class ConformSeries(SimpleInterface):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since the series is being mapped over, perhaps "ConformImage" or possibly "Conform"?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conform sounds good!

@@ -279,6 +228,93 @@ def _run_interface(self, runtime):
return runtime


class ConformSeriesInputSpec(BaseInterfaceInputSpec):
t1w = File(exists=True, mandatory=True,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reasons for using t1w_list are now obviated, so I'd switch to in_file, to make this more generic.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will do

target_shape = traits.List(traits.Int,
desc='Target shape information')
target_span = traits.List(traits.Float,
desc='Target span information')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, I'd drop target_span.



class ConformSeriesOutputSpec(TraitedSpec):
t1w = File(exists=True, desc='conformed T1w image')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

out_file = File(...)

@effigies
Copy link
Member

@anibalsolon No rush, but let us know if there's anything you need from us.

@anibalsolon
Copy link
Contributor Author

anibalsolon commented Aug 19, 2017

Hi @effigies, sorry for the delay!

Thank you for the detailed review. If there is something I can improve, please let me know.

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few quick comments.

Also, have you tested this on any of your own datasets with multiple T1w images?

while valid.any():
target_zooms = all_zooms[valid].min(axis=0)
scales = all_zooms[valid] / target_zooms
if np.all(scales < max_scale):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since it's only used once, don't bother setting max_scale = self.inputs.max_scale above. Just use self.inputs.max_scale here.

desc='input T1w images')

max_scale = traits.Float(3.0, usedefault=True,
desc='Maximum scaling factor in images to accept')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these can each be on a single line (max line length: 99), but I haven't checked.

In any event, we don't need a blank space between t1w_list and max_scale.

These comments apply to all of the Input/OutputSpecs.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice! Do you know if it is possible to enforce this pattern of attr spacing via flake?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not that I know of.

out_report = File(exists=True, desc='conformation report')


CONFORMATION_TEMPLATE = """\t\t<h3 class="elem-title">Anatomical Conformation</h3>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's keep these above the input/output specs. So the pattern is:

CONSTANTS

class XInputSpec(...):
    ...

class XOutputSpec(...):
    ...

class X(...):
    ...

#. Orient to RAS (left-right, posterior-anterior, inferior-superior)
#. Along each dimension, resample to minimum voxel size, maximum number of voxels
class TemplateDimensions(SimpleInterface):
"""Filter a series of T1w images based on the requirements for up-sampling.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Proposed replacement for this first line, focusing on the main goal of this interface:

Finds template target dimensions for a series of T1w images, filtering low-resolution images, if necessary.

Along each axis, the minimum voxel size (zoom) and the maximum number of voxels (shape) are found
across images.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great 👌

else:
copyfile(fname, out_name, copy=True, use_hardlink=True)

self._results['out_file'] = out_name
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's not bother copying, just point to the original file:

def _run_interface(self, runtime):
    # Load image, orient as RAS
    fname = self.inputs.in_file
    orig_img = nb.load(fname)
    reoriented = nb.as_closest_canonical(orig_img)

    if reoriented is not orig_img:
        out_name = fname_presuffix(fname, suffix='_ras', newpath=runtime.cwd)
        reoriented.to_filename(out_name)
    else:
        out_name = fname

    self._results['out_file'] = out_name

if reoriented is not orig_img:
reoriented.to_filename(out_name)
else:
copyfile(fname, out_name, copy=True, use_hardlink=True)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same thing here, wrt copying.

@anibalsolon
Copy link
Contributor Author

I tested on Poldrack's MyConnectome dataset, as suggested in issue #613
http://myconnectome.org/wp/data-sharing/ Sessions 13-24

Copy link
Member

@effigies effigies left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great! Thanks. One last minor comment on the docs, but I think this is good to go. Could you go ahead and merge master? I'll run it on a couple subjects we've had issues with in the past to make sure there are no regressions.

Performs two basic functions:

#. Orient to RAS (left-right, posterior-anterior, inferior-superior)
#. Along each dimension, resample to minimum voxel size, maximum number of voxels
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is no longer calculating the zooms, shape, we can just say:

"resample to target zooms (voxel sizes) and shape (number of voxels)"

@anibalsolon
Copy link
Contributor Author

Awesome!

I cannot merge because I don't have write access.

If there are some datasets that I can help testing, just let me know.

@chrisgorgo
Copy link
Contributor

I think @effigies meant merging "master" branch from the upstream (poldracklab/fmriprep) into your branch.

@anibalsolon
Copy link
Contributor Author

Ok! Merged with the upstream. Thanks @chrisfilo

@effigies
Copy link
Member

Good to go. Thanks a lot!

@effigies effigies merged commit 9542c4f into nipreps:master Aug 21, 2017
@effigies effigies mentioned this pull request Sep 20, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants