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

Affine registration PR 3/3 #654

Merged
merged 47 commits into from Jul 26, 2015

Conversation

Projects
None yet
7 participants
@omarocegueda
Contributor

omarocegueda commented May 24, 2015

This is the last PR of the affine registration series. It also includes a tutorial explaining how to do an ANTS2-like affine registration. The following graph shows the performance of this implementation compared to ANTS2, using sparse sampling (30% of the voxels).
image

EDIT: we have now the same accuracy as ANTS2 in this dataset (image updated). Our implementation is a bit faster too (still a bit slower than ANTS2), the following are the updated times:

(sparse means use 30% of the voxels for PDF estimation)
ANTS2 sparse: 11.0 minutes
Dipy sparse: 12.0 minutes
Flirt : 4.41 minutes
Nipy: 49 seconds

@@ -96,310 +97,6 @@ def get_direction_and_spacings(affine, dim):
A = affine[:dim,:dim]
return A.dot(np.diag(1.0/scalings)), scalings
class ScaleSpace(object):

This comment has been minimized.

@omarocegueda

omarocegueda May 24, 2015

Contributor

I moved the ScaleSpace class to its own module because now we have two different ways of building a scale space (the ANTS-like anisotropic scale space and the ANTS2-like isotropic one). The new scale space is IsotropicScaleSpace.

@@ -1219,10 +1219,10 @@ cdef double _compute_mattes_mi(double[:, :] joint,
if mi_gradient is not None:
for k in range(n):
mi_gradient[k] -= joint_gradient[i, j, k] * factor

This comment has been minimized.

@omarocegueda

omarocegueda May 24, 2015

Contributor

In the original paper by Mattes, he uses the negative MI, and I coded everything like in the paper, but it may be misleading to the users if we return the negative MI. The change of sign (to do minimization instead of maximization) will be made at the last moment, after we call this from the optimizer.

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

That's the -1 * on line 258 of immafine.py? Just making sure that I understand this. Is it worth adding comment here, so that we don't forget?

import scipy.ndimage.filters as filters
class ScaleSpace(object):

This comment has been minimized.

@omarocegueda

omarocegueda May 24, 2015

Contributor

This class is to build an ANTS-like scale space. Here, the sigma along each x,y,z direction depends on the voxel size (this SS is anisotropic)

return self._get_attribute(self.sigmas, level)
class IsotropicScaleSpace(ScaleSpace):

This comment has been minimized.

@omarocegueda

omarocegueda May 24, 2015

Contributor

And this is to build an ANTS2-like SS. This provides the user with more flexibility to control how the SS is built, although it may be harder to use...

@@ -315,29 +315,6 @@ def test_optimizer_exceptions():
assert_raises(ValueError, optimizer._get_energy_derivative)
def test_scale_space_exceptions():

This comment has been minimized.

@omarocegueda

omarocegueda May 24, 2015

Contributor

The test was moved to test_scalespace.py.

"""
Finally, lets refine with a full affine transform (translation, rotation, scale and
shear), it is safer to fit more degrees of freadom now, since we must be very close

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

"freadom" => "freedom"

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

fixed =)

opt_tol=opt_tol,
factors = factors,
sigmas = sigmas)

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Inconsistent use of spaces

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

fixed

histograms
sampling_proportion : int in (0,100]
the percentage of voxels to be used for estimating the (joint and
marginal) intensity histograms. If None, dense sampling is used.

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Default value is 100? I think that "dense sampling" is not specific enough.

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

You are right!, I forgot to mention. I added the following note:

        Notes
        -----
        Since we use linear interpolation, images are not, in general,
        differentiable at exact voxel coordinates, but they are differentiable
        between voxel coordinates. When using sparse sampling, selected voxels
        are slightly moved by adding a small random displacement within one
        voxel to prevent sampling points to be located exactly at voxel
        coordinates. When using dense sampling, this random displacement is
        not applied.
get_direction_and_spacings(moving_grid2space, self.dim)
self.prealign = prealign
self.param_scales = None
static_32 = np.array(static).astype(np.float32)

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

These _32 arrays only get used in the case where sampling_proportion is not none. Maybe move them into that else clause?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Of course! thanks! =)

the parameter vector of the transform currently used by the metric
(the transform name is provided when self.setup is called), n is
the number of parameters of the transform
update_gradient : Boolean

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Optional + what is the default value.

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Good to know! I never specified optional previously, now I see what PEP8 suggests

This comment has been minimized.

@arokem

arokem May 26, 2015

Member

Best place to look is in the numpy docstring standard:

https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt

On Mon, May 25, 2015 at 7:52 PM, Omar Ocegueda notifications@github.com
wrote:

In dipy/align/imaffine.py
#654 (comment):

  •    r""" Updates marginal and joint distributions and the joint gradient
    
  •    The distributions are updated according to the static and warped
    
  •    images. The warped image is precisely the moving image after
    
  •    transforming it by the transform defined by the xopt parameters.
    
  •    The gradient of the joint PDF is computed only if update_gradient
    
  •    is True.
    
  •    Parameters
    

  •    xopt : array, shape (n,)
    
  •        the parameter vector of the transform currently used by the metric
    
  •        (the transform name is provided when self.setup is called), n is
    
  •        the number of parameters of the transform
    
  •    update_gradient : Boolean
    

Good to know! I never specified optional previously, now I see what PEP8
suggests


Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/654/files#r31001768.

moving_32 = np.array(moving).astype(np.float32)
T = np.eye(self.dim + 1)
if self.prealign is not None:

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Once the prealignment happens, do you need to set it to None? Or to np.eye? Will it subseqeuntly get reapplied otherwise (e.g. below in update_dense)?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

There's no risk of reapplying the prealign matrix because images are never warped more than once: each time we need to do warping, the full composition is recomputed and applied to the original image.

an instance of a metric
level_iters : list
the number of iterations at each level of the Gaussian pyramid.
level_iters[0] corresponds to the finest level, level_iters[n-1]

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Add single quotes around level_iters[0] and level_iters[n-1]. Also, that can probably be level_iters[-1], no?

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Also, that's not what the example says: there you say that the last item in this list corresponds to the finest resolution. Here you say it corresponds to the coarsest. So, which is it?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Oh! that's right!, it's a "legacy documentation": at the beginning 0-th level was the finest, but then we inverted the order to make it consistent with the order in which optimization is performed (coarsest first, finest last), nice catch!

method : string
optimization method to be used
ss_sigma_factor : float
parameter of the scale-space smoothing kernel. For example, the

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

This is a 3D Gaussian?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Yes (or 2D for 2D images). I added Gaussian to the docstring

parameter of the scale-space smoothing kernel. For example, the
std. dev. of the kernel will be factor*(2^i) in the isotropic case
where i = 0, 1, ..., n_scales is the scale
options : None or dict,

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Why is opt_tol a separate parameter from the options?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Yeah, sorry, this is ugly, I am not even using it! =(. I think it's better to put it in options.

The long explanation is that, at the beginning, I was not using the Optimizer class, but I implemented my own preconditioned conjugate gradient with golden section line search, and then I thought that nobody was going to accept that, given that we already have the optimization tools, so I removed it and used the Optimizer class instead.

factors : list of floats
custom scale factors to build the scale space (one factor for each
scale)
method : string

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Same question here as below: why are some things part of the options dict and some things explicitly specified?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

In this case, method is an explicit argument of the Optimizer class:
https://github.com/nipy/dipy/blob/master/dipy/core/optimize.py#L24
even though it receives options too

Parameters
----------
xopt : array, shape (n,)

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Returns section in the docstring

----------
xopt : array, shape (n,)
the parameter vector of the transform currently used by the metric
(the transform name is provided when self.setup is called), n is

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Returns section in this docstring

elif prealign == 'centers':
self.prealign = aff_geometric_centers(static, static_grid2space,
moving, moving_grid2space)
else:

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Do you want to check whether this is an array, before setting it? If someone typos one of the other options (e.g. types masss instead of mass), you probably want to raise an error, instead of silently assigning this string to self.prealign

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Right! I added
elif isinstance(prealign, np.ndarray) and prealign.shape == (n,):
to check for number of parameters passed too. Else, raise ValueError.

static image: it will provide the grid and grid-to-space transform for
the warped image
static_grid2space:
grid-to-space transform associated to the static image

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

=> "associated with"

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Oh! no! again! what's wrong with me!? please don't tell Matthew! =S

moving: array, shape(S', R', C')
moving image
moving_grid2space:
grid-to-space transform associated to the moving image

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Same here:

=> "associated with"

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Yeah, sorry, I just fixed many more in this file and scalespace too

Returns
-------
transform : array, shape(4, 4)
the affine transformation (translation only, in this case) aligning

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

What happens when the two images are sampled at different resolutions? Is this still just a translation?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Yes, all transformations are represented in physical space (the actual transform mapping voxels to voxels is computed only when warping is done), so aligning the centers of mass is always a translation regardless of the images' sampling.

Gaussian kernel with increasing smoothing parameter. If the image's
voxels are isotropic, the smoothing will be the same along all
directions: at level L = 0,1,..., the sigma is given by
s * ( 2^L - 1 ).

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Maybe use a math sphinx directive here, or mathdollars

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Done

weaker along low resolution directions.
Parameters
----------

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Please specify which parameters are optional and what their default values are.

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Done, will do a full check to make sure I specify all optionals

# Compute the rest of the levels
min_spacing = np.min(input_spacing)
for i in range(1, num_levels):
scaling_factor = 2**i

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

White-space around the operator

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Fixed

scale and smoothing factors for all scales.
Parameters
----------

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

Again, which are optional and what are their defaults (especially where the kwarg default value is None, it's not obvious what that would entail!)

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Done

[0, 0, 0, 1]])
trans_inv = np.linalg.inv(trans)
for theta in [-1 * np.pi/6.0, 0.0, np.pi/5.0]: #rotation angle

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

What's wrong with the following?

 for rotation_angle in [-1 * np.pi/6.0, 0.0, np.pi/5.0]:
      for scale_factor in [0.83,  1.3, 2.07]:
             ...

That way, you don't need these comments, because the variable name explains what it is

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

Yeah, the main reason was to avoid breaking the lines below (because of the line length limit) and to make the below matrix look more compact. That's the PEP8 rule I like the least: it makes it harder to choose variable names, especially in cython with the nested loops. Now using rotation_angle and scale_factor.

start_sad = np.abs(static - moving).sum()
metric = imaffine.MattesMIMetric(32, sampling_pc)
affreg = imaffine.AffineRegistration(metric,
[10000, 1000, 100], 1e-5,

This comment has been minimized.

@arokem

arokem May 25, 2015

Member

These are all default settings, right? What happens when you use other settings?

This comment has been minimized.

@omarocegueda

omarocegueda May 26, 2015

Contributor

These are just settings to obtain a "good" registration so we can verify that the final alignment is good. Let me think what other settings would be interesting to test ("malicious" settings like lists of strings and things like that, maybe?).

@arokem

This comment has been minimized.

Member

arokem commented May 25, 2015

Seems to currently be failing on the minimal scipy requirement (0.9), because the setting of inputs to the optimization have changed since:

https://travis-ci.org/nipy/dipy/jobs/63966654#L2007

Any way we can provide some kind of compatibility shim in dipy.core.optimize to help with that?

@arokem

This comment has been minimized.

Member

arokem commented May 25, 2015

Any chance to increase test coverage in imaffine? It's currently at 90% (which is great), but I suspect that some non-default settings are not getting exercised at all.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented May 30, 2015

Hi @Garyfallidis, @arokem, @matthew-brett,
I think this is a bug in dipy.core.optimize for scipy<=0.12:
https://github.com/nipy/dipy/blob/master/dipy/core/optimize.py#L169
According to this:
http://docs.scipy.org/doc/scipy-0.9.0/reference/generated/scipy.optimize.fmin_l_bfgs_b.html#scipy.optimize.fmin_l_bfgs_b
this function expects fprime as third argument, but is given the tuple args instead, which is not callable and causes this error when it attempts to get the gradient:
https://travis-ci.org/nipy/dipy/jobs/63966654#L2021
As a quick fix, I am simply passing jac:
omarocegueda@6dbcd4a#diff-6b28d8cc8891944f1a6adca73299277aR154
but I think this may need a bit more thought about what's a good way to keep compatibility with scipy<=0.12. What do you think?
Thanks! =)

@arokem

This comment has been minimized.

Member

arokem commented Jul 27, 2015

Oh durnit - I relaunched that run, and now it's working.

It was failing on this assertion:

https://github.com/nipy/dipy/blob/master/dipy/align/tests/test_imaffine.py#L307

But only on python 3.3 (?!). Let's hope that was just a fluke.

On Mon, Jul 27, 2015 at 12:45 PM, Omar Ocegueda notifications@github.com
wrote:

Hi @arokem https://github.com/arokem,
I don't see the failure, is that the right link?


Reply to this email directly or view it on GitHub
#654 (comment).

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 27, 2015

Thanks for reporting that! this kind of weird failures usually indicates a memory access violation, or an uninitialized buffer. I'll double check.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 27, 2015

No need to revert, we're just discussing the details of the inheritance, and I doubt anyone will start using that stuff in the next week or so.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Hi @arokem,

Oh durnit - I relaunched that run, and now it's working.

Marc just reported the same failure in Python 3.4. I am investigating but since this is already in master, do you think it is preferable to skip the test for now and open an issue so it stops breaking the tests (I can re-activate the test after I figure out what the issue is about)?

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

Yeah - I saw that. I am somewhat reluctant to remove this from our tests,
because it's better to have a constant reminder of this, than to forget
about it in the shuffle of doing other things... Do you have a clear path
forward to fix this? Any idea why this didn't show up during the review on
this PR? I haven't looked closely at the code that leads to this yet.

On Wed, Jul 29, 2015 at 7:21 AM, Omar Ocegueda notifications@github.com
wrote:

Hi @arokem https://github.com/arokem,

Oh durnit - I relaunched that run, and now it's working.

Marc just reported the same failure in Python 3.4. I am investigating but
since this is already in master, do you think it is preferable to skip the
test for now and open an issue so it stops breaking the tests (I can
re-activate the test after I figure out what the issue is about)?


Reply to this email directly or view it on GitHub
#654 (comment).

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

It's hard to reproduce (that's what I am trying to do right now), I think that's why it didn't show up during the review, it works correctly in my Ubuntu, Windows 8.1 and Travis didn't complain either.

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

The worst kind... I will also give it a try.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

I see this now too:

http://nipy.bic.berkeley.edu/builders/dipy-py2.6-32/builds/555/steps/shell_6/logs/stdio

dipy.align.tests.test_imaffine.test_align_origins_3d ... ok
dipy.align.tests.test_imaffine.test_affreg_all_transforms ... 
command timed out: 1200 seconds without output, attempting to kill
process killed by signal 9

Maybe that would be a good machine to start on?

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

@matthew-brett - is that really the same issue? Superficially it looks like another test than this failure. Not saying we don't need to fix this one as well - just trying to understand.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

Actually, even better, it's failing in a similar way on OSX, with a segfault:

http://nipy.bic.berkeley.edu/builders/dipy-py2.7-osx-10.8/builds/377/steps/shell_6/logs/stdio

That machine is much faster.

Same on my laptop:

$ nosetests dipy/align/tests/test_imaffine.py:test_affreg_all_transforms
nose.config: INFO: Ignoring files matching ['^\\.', '^_', '^setup\\.py$']
dipy.align.tests.test_imaffine.test_affreg_all_transforms ... Segmentation fault: 11
@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

Ariel - I'm guessing there's some memory access violation somewhere, maybe causing both problems, just recommending this one as a place to start because it's easier to reproduce, and hoping that this will lead us to a fix for both issues.

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

Actually - I get the test_mi_gradient failure on my laptop as well. At least for one run out one. Now trying to see if it's an intermittent thing.

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

*one run out of one

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Oh! that's awesome!, @matthew-brett, when you say "Maybe that would be a good machine to start on" do you mean using try-branch for debugging or is there a way to get access to the actual machine?

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

Sadly, intermittent: one for two.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

I can't give you access to the '10.8' machine (actually, it's a 10.10 machine now). But I can give you access to the 10.7 machine that is failing in the same way. I'm afraid that machine is in really bad shape, running very slowly. You should now have access via your ssh-rsa key for jomaroceguedag@gmail.com : ssh buildslave@169.229.158.3 - if that gets too painful, we can ask Min, who owns the faster machine, if he will give you access.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

That machine already had the code compiled. To get going using the virtualenv that the tests used:

$ cd osx-10.7/dipy-py2_7-osx-10_7/build/
$ source venv/bin/activate
@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Thanks @matthew-brett!, I just activated the venv and I am now trying to build with make ext from
/Users/buildslave/osx-10.7/dipy-py2_7-osx-10_7/source
But gcc fails with:
lipo: can't figure out the architecture type of: /var/folders/8n/t5rvqnld23n0mn_lbhkgyg4w0000gt/T//ccvZTm15.out error: command 'gcc-4.2' failed with exit status 1

is that the right thing to do?

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

You can get the parameters that buildbot runs with from looking at the file /Library/LaunchDaemons/edu.berkeley.bic.pimba.osx-10.7.plist. In this case I think you need:

export CC=clang
export MACOSX_DEPLOYMENT_TARGET=10.6
make ext
@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Awesome! it's compiling now =)

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Hi @matthew-brett, @arokem,
The segmentation fault is caused by the interpolation functions, e.g.:
https://github.com/nipy/dipy/blob/master/dipy/align/vector_fields.pyx#L201
Here, I am assuming that the coordinates are valid floating point values (not nan). Infinite is handled correctly. The problem is that at some point the optimizer requests the metric to be evaluated with a ScaleTransform with parameter nan, which creates an affine containing nans. When we ask numpy to invert this matrix with np.linalg.inv, numpy returns a matrix full of nans instead of raising an exception. Since the exception is not caught, I assume that it's a valid transform and attempt to transform the moving image with the provided matrix, which calls the interpolation function with parameter nan.

We have two options here:

  1. Add comparisons of the form "if x!=x return" to all interpolation functions , which will cause some overhead on all the registration module because these functions are called once per voxel each time we interpolate. I will need to actually measure how large is this overhead.
  2. Check that the matrix does not contain nans at the beginning of all the warping functions. The overhead will likely be much lower than (1).
  3. Raise an exception from imaffine when the requested parameter is invalid. This is a very simple, 2-line fix, but we will still assume that the input arguments to the module arevalid, and will probably get the "segfault" again if the user provides an invalid matrix, or displacement field.

I personally would vote (2). What do you think?

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

On Wed, Jul 29, 2015 at 11:00 AM, Omar Ocegueda notifications@github.com
wrote:

Hi @matthew-brett https://github.com/matthew-brett, @arokem
https://github.com/arokem,
The segmentation fault is caused by the interpolation functions, e.g.:
https://github.com/nipy/dipy/blob/master/dipy/align/vector_fields.pyx#L201
Here, I am assuming that the coordinates are valid floating point values
(not nan). Infinite is handled correctly. The problem is that at some point
the optimizer requests the metric to be evaluated with a ScaleTransform
with parameter nan, which creates an affine containing nans. When we ask
numpy to invert this matrix with np.linalg.inv, numpy returns a matrix full
of nans instead of raising an exception. Since the exception is not caught,
I assume that it's a valid transform and attempt to transform the moving
image with the provided matrix, which calls the interpolation function with
parameter nan.

We have two options here:

  1. Add comparisons of the form "if x!=x return" to all interpolation
    functions , which will cause some overhead on all the registration module
    because these functions are called once per voxel each time we interpolate.
    I will need to actually measure how large is this overhead.
  2. Check that the matrix does not contain nans at the beginning of all the
    warping functions. The overhead will likely be much lower than (1).
  3. Raise an exception from imaffine when the requested parameter is
    invalid. This is a very simple, 2-line fix, but we will still assume that
    the input arguments to the module arevalid, and will probably get the
    "segfault" again if the user provides an invalid matrix, or displacement
    field.

I personally would vote (2). What do you think?

Yes, unless I am missing something, 2 definitely sounds like the best way
to go about this.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jul 29, 2015

How about doing 2 and 3? Is there any reason not to raise the error for option 3?

@arokem

This comment has been minimized.

Member

arokem commented Jul 29, 2015

P.S. Does this explain why the error is intermittent? With the random seed
fixed at the beginning of this test, that's not obvious to me.

On Wed, Jul 29, 2015 at 11:09 AM, Ariel Rokem arokem@gmail.com wrote:

On Wed, Jul 29, 2015 at 11:00 AM, Omar Ocegueda notifications@github.com
wrote:

Hi @matthew-brett https://github.com/matthew-brett, @arokem
https://github.com/arokem,
The segmentation fault is caused by the interpolation functions, e.g.:
https://github.com/nipy/dipy/blob/master/dipy/align/vector_fields.pyx#L201
Here, I am assuming that the coordinates are valid floating point values
(not nan). Infinite is handled correctly. The problem is that at some point
the optimizer requests the metric to be evaluated with a ScaleTransform
with parameter nan, which creates an affine containing nans. When we ask
numpy to invert this matrix with np.linalg.inv, numpy returns a matrix full
of nans instead of raising an exception. Since the exception is not caught,
I assume that it's a valid transform and attempt to transform the moving
image with the provided matrix, which calls the interpolation function with
parameter nan.

We have two options here:

  1. Add comparisons of the form "if x!=x return" to all interpolation
    functions , which will cause some overhead on all the registration module
    because these functions are called once per voxel each time we interpolate.
    I will need to actually measure how large is this overhead.
  2. Check that the matrix does not contain nans at the beginning of all
    the warping functions. The overhead will likely be much lower than (1).
  3. Raise an exception from imaffine when the requested parameter is
    invalid. This is a very simple, 2-line fix, but we will still assume that
    the input arguments to the module arevalid, and will probably get the
    "segfault" again if the user provides an invalid matrix, or displacement
    field.

I personally would vote (2). What do you think?

Yes, unless I am missing something, 2 definitely sounds like the best way
to go about this.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

Oh! of course!, we need 2 and 3 so we can catch the error and send inf to the optimizer when it attempts to evaluate the metric at an invalid point. We don't want to raise the error because we want the optimizer to continue exploring even when it evaluates the metric at an invalid point.

@arokem, I'm not sure if this explains the intermitent behavior. I think it depends on what happens when we attempt to access an array at index 'nan'. Let me investigate a bit more about that, but for now this should fix this memory error too:
http://nipy.bic.berkeley.edu/builders/dipy-bdist32-33/builds/138/steps/shell_8/logs/stdio

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 29, 2015

what happens when we attempt to access an array at index 'nan'

Sorry, I meant: access an array at the index that results from converting nan to integer (the resulting index might yield a valid memory address some times and some times it might not)...

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Aug 3, 2015

Hello @omarocegueda , I noticed that you haven't renamed the following functions in what is now in the master
In [1]: from dipy.align.imaffine import align_
align_centers_of_mass align_geometric_centers align_origins

Can you please use transform_X instead as we discussed in this PR?
Thx in advance.

@coveralls

This comment has been minimized.

coveralls commented Apr 5, 2017

Coverage Status

Changes Unknown when pulling 741e3c0 on omarocegueda:imaffine_pr3 into ** on nipy:master**.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment