# Optimal transport warping implementation in Matlab

I am implementing the paper "Optimal Mass Transport for Registration and Warping", my goal being to put it online as I just cannot find any eulerian mass transportation code online and this would be interesting at least for the research community in image processing.

The paper can be summarized as follows :
- find an initial map $u$ using 1D histogram matchings along the x and y coordinates
- solve for the fixed point of $u_t = \frac{1}{\mu_0} Du \nabla^\perp\triangle^{-1}div(u^\perp)$ , where $u^\perp$ stands for a 90 degrees counter clockwise rotation, $\triangle^{-1}$ for the solution of the poisson equation with Dirichlet boundary conditions (=0), and $Du$ is the determinant of the Jacobian matrix.
- stability is guaranteed for a timestep $dt<\min|\frac{1}{\mu_0}\nabla^\perp\triangle^{-1}div(u^\perp)|$

For numerical simulations (performed on a regular grid), they indicate using matlab's poicalc for solving the poisson equation, they use centered finite differences for spatial derivatives, except for $Du$ which is computed using an upwind scheme.

Using my code, the energy functional and curl of the mapping are properly decreasing for a couple iterations (from a few tens to a few thousands depending on the time step). But after that, the simulation explodes : the energy increases to reach a NAN in very few iterations. I tried several orders for the differentiations and integrations, and different interpolation schemes, but I always get the same issue (even on very smooth images, non-zero everywhere etc.).
Anyone would be interested in looking at the code and/or the theoretical problem I am facing ? The code is rather short.

I am only interested in the optimal transport part of the paper for now, not the additional regularization term.

Thanks !

edit: a higher order integration code to replace the cumptrapz can be found here
edit2: please replace the call to gradient2() at the end, to gradient(). This was actually a higher order gradient but doesn't solve things.

-
Though I do not believe this question is completely off-topic here, you may have more success trying to ask it at scicomp.stackexchange.com – András Bátkai May 8 '12 at 5:32
Thanks for the tip - I didn't know this one and was indeed looking for something between mathoverflow and stackoverflow, with a higher emphasis on math than computer science. If I don't get relevant answers here, I'll post it on scicomp. Thanks ! – WhitAngl May 8 '12 at 6:28
You are aware of the "earth movers' distance", I presume? – Felix Goldberg May 8 '12 at 8:30
yes I am :) But all EMD codes online are using the linear programming formulation rather than the PDE. – WhitAngl May 8 '12 at 14:12
as suggested by András, I forwarded the question to scicomp.stackexchange : scicomp.stackexchange.com/questions/2162/… . Thks. – WhitAngl May 9 '12 at 18:59