2
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ 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 $\endgroup$ Commented May 8, 2012 at 5:32
  • 1
    $\begingroup$ 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 ! $\endgroup$ Commented May 8, 2012 at 6:28
  • $\begingroup$ You are aware of the "earth movers' distance", I presume? $\endgroup$ Commented May 8, 2012 at 8:30
  • $\begingroup$ yes I am :) But all EMD codes online are using the linear programming formulation rather than the PDE. $\endgroup$ Commented May 8, 2012 at 14:12
  • $\begingroup$ as suggested by András, I forwarded the question to scicomp.stackexchange : scicomp.stackexchange.com/questions/2162/… . Thks. $\endgroup$ Commented May 9, 2012 at 18:59

0

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.