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 2/3 #628

Merged
merged 18 commits into from May 24, 2015

Conversation

Projects
None yet
4 participants
@omarocegueda
Contributor

omarocegueda commented Apr 13, 2015

This is the second PR to add affine registration to the align module. It contains the low level functions to compute the Mattes' Mutual Information metric and its analytic gradient w.r.t. the parameters of the transformations defined in PR1.

class MattesBase(object):
def __init__(self, nbins):
r""" MattesBase
Base class for the Mattes' Mutual Information metric.

This comment has been minimized.

@arokem

arokem Apr 17, 2015

Member

Reference?

This comment has been minimized.

@arokem

arokem Apr 17, 2015

Member

Also: Mattes's, if this is a possessive, or Mattes (no apostrophe) if not. I think the latter.

This comment has been minimized.

@omarocegueda

omarocegueda Apr 19, 2015

Contributor

Fixed =) (removed the apostrophe and added missing reference)

Just for the 'trivia', (don't put much attention to this, I have no preference about it =P ), apparently there's no rule for using apostrophe to indicate possession when the proper noun ends with 's', they just recommend to stay consistent:
http://www.grammarbook.com/punctuation/apostro.asp

This comment has been minimized.

@arokem

arokem Apr 19, 2015

Member

I go with Strunk and White: http://www.bartleby.com/141/strunk.html

On Sun, Apr 19, 2015 at 9:22 AM, Omar Ocegueda notifications@github.com
wrote:

In dipy/align/mattes.pyx
#628 (comment):

  •                                  _apply_affine_3d_x1,
    
  •                                  _apply_affine_3d_x2,
    
  •                                  _apply_affine_2d_x0,
    
  •                                  _apply_affine_2d_x1)
    
    +from dipy.align.transforms cimport (Transform)
    +
    +cdef extern from "dpy_math.h" nogil:
  • double cos(double)
  • double sin(double)
  • double log(double)

+class MattesBase(object):

  • def init(self, nbins):
  •    r""" MattesBase
    
  •    Base class for the Mattes' Mutual Information metric.
    

Fixed =) (removed the apostrophe and added missing reference)

Just for the 'trivia', (don't put much attention to this, I have no
preference about it =P ), apparently there's no rule for using apostrophe
to indicate possession when the proper noun ends with 's', they just
recommend to stay consistent:
http://www.grammarbook.com/punctuation/apostro.asp


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

# intensities. Padding is the number of extra bins used at each side
# of the histogram (a total of [2 * padding] extra bins). Since the
# support of the cubic spline is 5 bins (the center plus 2 bins at each
# side) we need a padding of 2, in this case.

This comment has been minimized.

@arokem

arokem Apr 17, 2015

Member

"this case" being the case where nbins is 1?

This comment has been minimized.

@omarocegueda

omarocegueda Apr 19, 2015

Contributor

sorry, I meant "the cubic spline case", I clarified it in the docstring.

Parameters
----------
nbins : int
the number of bins of the joint and margianl probability density

This comment has been minimized.

@arokem

arokem Apr 17, 2015

Member

typo: margianl => marginal

This comment has been minimized.

@omarocegueda

omarocegueda Apr 19, 2015

Contributor

Fixed

----------
nbins : int
the number of bins of the joint and margianl probability density
functions (the actual number of bins of the joint PDF is nbins^2)

This comment has been minimized.

@arokem

arokem Apr 17, 2015

Member

nbins^2 => $nbins^2$ OR nbins**2

This comment has been minimized.

@omarocegueda

omarocegueda Apr 19, 2015

Contributor

Fixed, used nbins**2. Sorry, this may be a silly question, but is the latex notation actually useful in the docstring, I mean, are the docstrings being parsed at some point to generate the math expressions?

This comment has been minimized.

@arokem

arokem Apr 19, 2015

Member

Yes - the API docs in the html rendering go through a latex rendering. For
example:
http://nipy.org/dipy/reference/dipy.reconst.html#fractional-anisotropy

On Sun, Apr 19, 2015 at 9:27 AM, Omar Ocegueda notifications@github.com
wrote:

In dipy/align/mattes.pyx
#628 (comment):

  •    appropriate optimizer.
    
  •    Notes: we need this class in cython to allow
    
  •    _joint_pdf_gradient_dense_2d and _joint_pdf_gradient_dense_3d to
    
  •    use a nogil Jacobian function (obtained from an instance of the
    
  •    Transform class), which allows us to evaluate Jacobians at all
    
  •    the sampling points (maybe the full grid) inside a nogil loop.
    
  •    The reason we need a class is to encapsulate all the parameters
    
  •    related to the joint and marginal distributions.
    
  •    Parameters
    

  •    nbins : int
    
  •        the number of bins of the joint and margianl probability density
    
  •        functions (the actual number of bins of the joint PDF is nbins^2)
    

Fixed, used nbins**2. Sorry, this may be a silly question, but is the
latex notation actually useful in the docstring, I mean, are the docstrings
being parsed at some point to generate the math expressions?


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

This comment has been minimized.

@omarocegueda

omarocegueda Apr 19, 2015

Contributor

It looks great! =D but is anything similar being done with the cython docs as well?

This comment has been minimized.

@arokem

arokem Apr 19, 2015

Member

Oh yeah - good point. I doubt it. At least I couldn't find any evidence of
it in our docs.

On Sun, Apr 19, 2015 at 12:05 PM, Omar Ocegueda notifications@github.com
wrote:

In dipy/align/mattes.pyx
#628 (comment):

  •    appropriate optimizer.
    
  •    Notes: we need this class in cython to allow
    
  •    _joint_pdf_gradient_dense_2d and _joint_pdf_gradient_dense_3d to
    
  •    use a nogil Jacobian function (obtained from an instance of the
    
  •    Transform class), which allows us to evaluate Jacobians at all
    
  •    the sampling points (maybe the full grid) inside a nogil loop.
    
  •    The reason we need a class is to encapsulate all the parameters
    
  •    related to the joint and marginal distributions.
    
  •    Parameters
    

  •    nbins : int
    
  •        the number of bins of the joint and margianl probability density
    
  •        functions (the actual number of bins of the joint PDF is nbins^2)
    

It looks great! =D but is anything similar being done with the cython docs
as well?


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Apr 17, 2015

@omarocegueda is this still WIP or we can go ahead and review?

@arokem

This comment has been minimized.

Member

arokem commented Apr 17, 2015

Oh, sorry, I missed that WIP marking and started commenting on some really
obvious stuff, while reading through and trying to wrap my head around
this. Feel free to ignore me :-)

On Friday, April 17, 2015, Eleftherios Garyfallidis <
notifications@github.com> wrote:

@omarocegueda https://github.com/omarocegueda is this still WIP or we
can go ahead and review?


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

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Apr 17, 2015

For me 'WIP' just means - don't merge yet - but do review.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Apr 17, 2015

I believe he said in the e-mail list that he should be ready around the 27th but I was not sure if that was for the 3rd PR or for this one.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Apr 17, 2015

Sorry guys!, I was offline and just saw your comments. I think this PR is ok to review now, I just would like to mention that there are 3 or 4 functions that are not tested yet, I am planning to add the missing tests this weekend. I think the 3rd PR should be much smaller because it will contain mostly python interfaces for the cython code of PRs 1 and 2. Thanks for your comments Ariel! =D.

@omarocegueda omarocegueda changed the title from WIP: affine registration PR 2/3 to Affine registration PR 2/3 Apr 19, 2015

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Apr 24, 2015

Hi @omarocegueda. Sorry for delaying to review this. I will do the review this Sunday.

moving = static.copy()
moving += np.random.randint(0, nvals//2, nr*nc).reshape(sh) - nvals//4
# This is just a simple way of making joint the distribution non-uniform

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

the joint

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

fixed =)

sh = (nr, nc)
static = np.random.randint(0, nvals, nr*nc).reshape(sh)
# This is just a simple way of making joint the distribution non-uniform

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

the joint

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

fixed

def create_random_image_pair_2d(nr, nc, nvals):
r""" Create a pair of images with an arbitrary, non-uniform joint

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

parameter description is missing

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Added missing documentation

# & Eubank, W. PET-CT image registration in the chest using
# free-form deformations. IEEE Transactions on Medical Imaging,
# 22(1), 120-8, 2003.
input = []

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

input is a reserved variable in Python. I would avoid using that.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Changed the name of the variable

assert_equal((indices%k).sum(), 0)
if __name__ == '__main__':

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

You don't need to list all the test functions here. You can just say np.testing.run_module_suite() and they will all run from nosetests.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Removed the function calls

transform parameters, and then communicate the results to the
appropriate optimizer.
Notes: we need this class in cython to allow

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

Notes should go at the end of the docstring and writen as in Parameters.

Notes
--------

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Moved notes to the bottom

class MattesBase(object):
def __init__(self, nbins):
r""" MattesBase
Base class for the Mattes Mutual Information metric [Mattes03].

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

Move this line up and remove MattesBase from the beginning of the string. It's redundant.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Done =)

cdef inline cnp.npy_intp _bin_index(double normalized, int nbins,
int padding) nogil:
r''' Index of the bin in which the normalized intensity 'normalized' lies

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member
``normalized`` lies.

This comment has been minimized.

@matthew-brett

matthew-brett Apr 29, 2015

Member

Need newline after first line for docstring standard.

I think you need single backticks around normalized to identify it as a parameter.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Good to know!, Using backticks now.

double sqrx = x * x
if absx < 1.0:
return ( 4.0 - 6.0 * sqrx + 3.0 * sqrx * absx ) / 6.0

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

No need for spaces next to parentheses

(4.0 - 6.0 * sqrx + 3.0 * sqrx *  absx) / 6.0

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Removed spaces

if absx < 1.0:
return ( 4.0 - 6.0 * sqrx + 3.0 * sqrx * absx ) / 6.0
elif absx < 2.0:
return ( 8.0 - 12 * absx + 6.0 * sqrx - sqrx * absx ) / 6.0

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

same here

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

Math correct 💃
Checked with referenced paper.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Great! thanks! =)

return -2.0 + 2.0 * x - 0.5 * x * x
else:
return 2.0 + 2.0 * x + 0.5 * x * x
return 0.0

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

Math checked! :)

mmarginal[j] = 0
for i in range(nbins):
mmarginal[j] += joint[i, j]

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 26, 2015

Member

pep 8 too many spaces

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Removed extra spaces. Sorry, I understand that it must be annoying to find so many PEP8 errors, so I did a full review using Spyder, as suggested.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Apr 26, 2015

I gave some initial comments. More reviewing tomorrow. Keep it up Omar!

self.joint_grad = None
self.metric_grad = None
self.metric_val = 0
self.joint = np.ndarray(shape = (self.nbins, self.nbins))

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

Clearer as np.zeros ?

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

Using np.zeros now =)

self.mmarginal = np.ndarray(shape = (self.nbins,), dtype = np.float64)
def bin_normalize_static(self, x):

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

Never called from Cython?

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

True, it's not being called from Cython, I just added the method for convenience in case the user needs it. It is being tested, though.

def update_pdfs_dense(self, static, moving, smask=None, mmask=None):
r''' Computes the Joint Probability Density Function of two images

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

Comment on where the JPDF goes? (Into self.smarginal, self.mmarginal I guess).

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

True, it's not obvious at all. The joint distribution goes to self.joint and the marginals to self.smarginal (marginal from the static image) and self.mmarginal (marginal from the moving image). Added comments explaining this.

raise ValueError(msg)
def update_pdfs_sparse(self, sval, mval):

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

By 'sparse' you mean, subsample pairs? Rename to reflect this? Like update_pdf_from_samples or something like that?

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 28, 2015

Member

I would document this better in the docstring to be more clear saying that this means subsampling etc. But the name of the function is quite nice. Dense vs sparse makes sense to me. I thought it was a good name when reviewing this part of the code.

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

That's right, sparse accepts a list of paired intensities sampled from both images at the same points in physical space. I extended the comments to make it more clear. Regarding the name, I have no strong preference, but if we change 'sparse' to 'from_samples' then we should also change 'dense' to 'from_grid' or something like that. I have a slight preference for 'sparse' and 'dense' because they're shorter. Do you think it is still confusing after extending the docstring?

self.mmarginal)
def update_gradient_dense(self, theta, transform, static, moving,

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

Sorry to read this quickly - but is this really an "update" - or does it do the whole calculation in one shot?

Say where the values go (?self.joint_grad)?

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

That's right, it does the full computation. Do you think something like "recompute" would be more clear?
Yes, it goes to self.joint_grad, I added a comment about that. =)

raise ValueError(msg)
def update_mi_metric(self, update_gradient=True):

This comment has been minimized.

@matthew-brett

matthew-brett Apr 27, 2015

Member

So this is an "update" if "update_gradient" is True, but not otherwise?

This comment has been minimized.

@omarocegueda

omarocegueda May 4, 2015

Contributor

It always computes the metric value, but the gradient will only be computed if "update_gradient" is True. This may save some time if the optimization algorithm doesn't need to compute the gradient but only its value at some point of the optimization (for example during line-search to determine the step size along the previously computed descent direction). I added a comment about it.

@omarocegueda omarocegueda force-pushed the omarocegueda:imaffine_pr2 branch from d529e08 to 422d8bb May 22, 2015

@arokem

This comment has been minimized.

Member

arokem commented May 22, 2015

If merging this into master on its own doesn't make sense, another option
is to make a PR against this branch, instead of against master. It might
make the review easier, because it will show only the diff relative to this
already completed work. And it will allow you to merge everything into
master together, so that we don't have a situation where somebody installs
from github and gets half-finished features.

On May 21, 2015 7:00 PM, "Eleftherios Garyfallidis" <
notifications@github.com> wrote:

Perfect!


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

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented May 22, 2015

That sounds good to me! (sorry, PR1 is already merged, so the unfinished
feature is already there, but it is a good exercise for future
split-features =) )

On Fri, May 22, 2015 at 10:53 AM, Ariel Rokem notifications@github.com
wrote:

If merging this into master on its own doesn't make sense, another option
is to make a PR against this branch, instead of against master. It might
make the review easier, because it will show only the diff relative to this
already completed work. And it will allow you to merge everything into
master together, so that we don't have a situation where somebody installs
from github and gets half-finished features.

On May 21, 2015 7:00 PM, "Eleftherios Garyfallidis" <
notifications@github.com> wrote:

Perfect!


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


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

"Cada quien es dueño de lo que calla y esclavo de lo que dice"
-Proverbio chino.
"We all are owners of what we keep silent and slaves of what we say"
-Chinese proverb.

http://www.cimat.mx/~omar

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented May 24, 2015

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 24, 2015

You don't like any of these options?

  1. X_grid2space X_space2grid
  2. X_img2space X_space2img
  3. X_img2scanner X_scanner2img
  4. X_affine and X_affine_inv
  5. X_grid2mm and X_mm2grid
  6. X_img2mm and X_mm2img
@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 24, 2015

So, 7. X_grid2world and X_world2grid right?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 24, 2015

In a recent tutorial that I wrote I used the same words
"grid to world" to explain the different coordinate systems. I think it's a good name too.
So, I am fine with 7 too.

    1. and 7. all good for me. Please make a decision and let's move to the next level.
@matthew-brett

This comment has been minimized.

Member

matthew-brett commented May 24, 2015

Yes, I prefer X_grid2world to the others.

I don't like the img2* options or X_affine

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented May 24, 2015

Alright!, thanks guys, I just changed the naming convention to X_grid2world, X_world2grid, and PR3+tutorial is ready. Should I wait until PR2 is merged or try to do PR3 on top of PR2, as Ariel suggested? what do you think?

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented May 24, 2015

if I do PR3 on top of PR2, I will be sending the PR to myself, right? (because PR2's branch is on my own fork)

Garyfallidis added a commit that referenced this pull request May 24, 2015

Merge pull request #628 from omarocegueda/imaffine_pr2
Affine registration PR 2/3

@Garyfallidis Garyfallidis merged commit f035833 into nipy:master May 24, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 24, 2015

Thx Omar. You can do PR3 on top of master.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented May 24, 2015

Awesome! thanks Elef!

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 24, 2015

GGT!

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented May 25, 2015

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