pyimreg / python-register

Linear and non-linear image registration methods, using scipy and numpy.
Other
62 stars 34 forks source link

Polynomial model #14

Open riaanvddool opened 13 years ago

riaanvddool commented 13 years ago

I would like to add a 2nd order polynomial model:

x' = a0 + a1.x + a2.y + a3.x.y + a4.x^2 + a5.y^2 y' = b0 + b1.x + b2.y + b3.x.y + b4.x^2 + b5.y^2

We use this quite a lot and it should be faster that the spline model because is only has 12 parameters.

nfaggian commented 13 years ago

Looks like this will be easy to derive the warp jacobian for - if you need help, let me know.

riaanvddool commented 13 years ago

Maybe you can check my attempt:

    n = self.coordinates.homogenous.shape[1]
    x = self.coordinates.homogenous[0,:]    
    y = self.coordinates.homogenous[1,:]

    dx = np.zeros((n, 12))
    dy = np.zeros((n, 12))

    dx[:,0] = 1.0
    dx[:,1] = x
    dx[:,2] = y
    dx[:,3] = x*y
    dx[:,4] = x**2
    dx[:,5] = y**2

    dy[:,6] = 1.0
    dy[:,7] = x
    dy[:,8] = y
    dy[:,9] = x*y
    dy[:,10] = x**2
    dy[:,11] = y**2

    return (dx, dy)
riaanvddool commented 13 years ago

I am getting the following error while trying out the Poly2 model:

Traceback (most recent call last): File "nonlinreg.py", line 60, in plotCB=plot.gridPlot File "/usr/local/lib/python2.6/dist-packages/python_register-0.1-py2.6-linux-i686.egg/register/register.py", line 225, in register p=p File "/usr/local/lib/python2.6/dist-packages/python_register-0.1-py2.6-linux-i686.egg/register/register.py", line 111, in __deltaP return np.dot( np.linalg.inv(H), np.dot(J.T, e)) File "/usr/lib/python2.6/dist-packages/numpy/linalg/linalg.py", line 355, in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype))) File "/usr/lib/python2.6/dist-packages/numpy/linalg/linalg.py", line 254, in solve raise LinAlgError, 'Singular matrix' numpy.linalg.linalg.LinAlgError: Singular matrix

nfaggian commented 13 years ago

I am having a look at this now - if the polynomial model is easy to invert it will be a great addition

nfaggian commented 13 years ago

Ok.

Riaan please have a look at how I have added the polynomial model you mentioned. I realized that there were a few problems.

1) in the optimization context identitiy should really map to no deformation - look at how the shift and affine warps add nothing to the warp. This parameterization was in Bakers paper.

2) The inverse warp is required for sampling - what I did in the spline model is a rough (and probably wrong approximation) :

 (grid coordinates minus warp) = inverse

We should clean up the examples before we merge this stuff - also I want to be sure that what I did is correct so please check the output of non-linreg.

While doing this I realized the comments are pretty out of date for some of the code - my fault :-(

riaanvddool commented 13 years ago

1) I though that with the polynomial model as i defined it, the parameters a1 and b2 should be =1 to get x= x and y = y?

2) Mm, not sure i understand.

nfaggian commented 13 years ago

1) For the optimization to work a starting point needs to be defined which is what I called "identity" in the code. The identity transform is the parameter set which yields no displacement in coordinates.

The deformation model (in all cases) is something like this:

x = x0 + dx y = y0 + dy

where x0, y0 are the grid coordinates and dx, dy are estimated using some sort of displacement model, beit affine or non-linear.

So if dx = f(p), the optimization relies on the fact that if p = 0, dx = 0.. same applies to dy.

2) For non-linear warps you have a nice equation for displacement (refer to the Kybic paper) but what is really computed (for sampling) is an inverse warp - otherwise the image will have holes in it after re-sampling. So I made an approximation (which is probably not great) and say that the inverse warp is:

x = x0 - dx y = y0 - dy

Hope this helps...

-N

riaanvddool commented 13 years ago

Please remember to try warping an image using the likeSpline model and look at the result...

nfaggian commented 13 years ago

Do we still want this model?

nfaggian commented 13 years ago

Made a branch for this work, ticket_0014