# Optimally "morph" one set of points into another

Given two sets of (two-dimensional) points, say, $$A=\{a_1,a_2,\ldots,a_{n_a}\}$$ and (you guessed it) $$B=\{b_1,b_2,\ldots,b_{n_b}\}$$, and $$d^2_{i,j}=\mid a_i-b_j\mid^{\ 2}$$ the "matrix" (not necessarily square) of distances between them, I want to map $$m:A\to B$$ (i.e., with each $$A$$-point $$a\in A$$ associate a corresponding $$B$$-point $$m(a)\in B$$) such that the following...

The maximal point-to-point $$d_{i,m(i)}$$ distance should be minimal. And once you've accomplished that first $$m:i\to m(i)$$ correspondence, each set gets one point smaller, $$A\to A\backslash\{a_i\}$$ and $$B\to B\backslash\{b_{m(i)}\}$$. And now, the second point-to-point correspondence made by $$m$$ should similarly minimize the maximal distance on these one-point-smaller sets. And then iteratively. So, my question: what's an efficient algorithm for this? ($$n_a,n_b\sim10^4$$, way too large for brute force)

And it's not necessary that $$n_a=n_b$$. If $$n_a\lt n_b$$, discard any of the excess $$n_b-n_a$$ $$B$$-points you like such that the set of all those maximal distances is minimized. That is, just do the iteration, and when you're done, just discard the "unused" $$B$$-points. On the other hand, if $$n_b\lt n_a$$ just discard the $$A$$-points resulting from the first $$n_a-n_b$$ iteration steps (i.e., again minimize the set of maximal distances). At least, I think that's the right "discard" scheme, but please correct me if I'm wrong.

The "point" (sorry:) of all this is, as per the subject, to "morph" one bitmapped image into another, whereby each $$A$$-point denotes a source image pixel, and each $$B$$-point denotes a target image pixel (and if $$n_b\gt n_a$$, I'll just inconspicuously/randomly "pop" a few extra $$B$$-pixels into place during each frame of the morph). Since this kind of morphing is all over the place, I'd have thought I could easily google algorithms/code/etc. But I strangely couldn't get google to cough it up. So if you're familiar with this kind of stuff and can just point me in the right direction, that would be great, too. Thanks.

E d i t
-----------------

Thanks again, @LedHead and @BrunoLevy. It's been quite some time (better late than never), but I finally got around to implementing this stuff, more or less as originally intended. Here's an illustration, followed by some explanatory remarks... The interpolation's nearest neighbor as suggested by you, one-to-many if the source image has fewer pixels than the target (and vice-versa if vice-versa). But the path from source-->target pixel isn't a straight line. Instead, it's a bezier curve whose first control point is the source row,col and whose last control point is the target's. In-between control points determine the shape of the path. They're basically chosen randomly, but constrained within the row,col-bounds of the image. However, when specifying the gif, you can choose multipliers for those constraints. For the illustrated gif, I chose constraints keeping all control points near the center. And that's why, in this case, the morphing first contracts all the pixels inwards, and then expands them outwards.

By the way, bezier also includes straight-line interpolation between nearest neighbors, simply by choosing just two control points, which are then necessarily source coordinates and target coordinates. The exact same gif, changed only in that way, looks like... • Does nearest neighbor interpolation not work for you? Your algorithm is more complicated than this, but I don't see the motivation for it. Jun 16, 2017 at 6:29
• Thanks @led23head I was certainly aware of "interpolation", but not "nearest neighbor interpolation". That wikipedia page doesn't seem to say very much, but google also coughed up mathworks.com/help/vision/ug/interpolation-methods.html that doesn't seem to suggest applicability to this problem. On the other hand, it's pretty intuitively obvious what you're suggesting, which is certainly applicable and may be adequate for my purposes (as described at end of comment to Bruno, below). And it's pretty simple to program, so I'll probbaly give it a try, just to see how it looks. Jun 16, 2017 at 22:30
• @led23head Oops, I'm seeing a slight complication: one target $B$-point may be the nearest neighbor of several source $A$-points. In that case I'm thinking to choose the $A$-point whose distance to it is largest, to avoid making the already-worst-case any "worser". Presumably, any very-nearby $A$-point can still find another pretty-near $B$-neighbor, i.e., nearer than the worst-case to which the original $B$-point was assigned. That'll be a little trickier to program, and possibly much more computationally intensive, but I'm not yet seeing how to work out all the details. Jun 16, 2017 at 22:55
• I think you want to think about it in reverse. For each point in $B$, find the nearest point in $A$. Is there some reason you can't have two $B$ points copy the same $A$ point? Jun 16, 2017 at 23:03
• Thanks again @led23head . Yeah, that may work better. But one-to-many (one_$A$-point$\to$many_$B$-points) may leave many of the $A$-points with "nowhere to go" (especially when $n_a\gt n_b$ to begin with). So the intermediate animation frames might not look like a "seamless morphing" to the viewer. Hard to tell without trying (which isn't entirely trivial), so I'm thinking, ab initio, the more one-to-one the mapping, the better. But it's ultimately probably a three-way-trade-off between programming effort, computational complexity, and viewer experience. I'd doubt any solution optimizes all. Jun 17, 2017 at 0:24

You may consider numerical optimal transport. It does not exactly fit your specification, but for your image morphing application it may be well suited. In the discrete setting, given your two set of points $A$ and $B$, optimal transport computes an assignment matrix $m$ between $A$ and $B$ that has unit row sums and unit column sums (a so-called bi-stochastic matrix) and that minimizes $\sum_i \sum_j m_{ij} c_{ij}$ for some costs $c_{ij}$ (for instance, $c_{ij}$ = the Euclidean distance between point $i$ and point $j$.

Optimal Transport is a general theory, that applies in differente settings (comprising continuous functions and discrete objects like your point sets). There are good books on the general theory  and how to use it in practice .

In your discrete case, the standard algorithm to be used is the so-called "auction algorithm" , but it is very slow. Thus, some accelerated counterparts were invented , based on the observation that adding a regularization term (entropy) makes it possible to use a much faster algorithm (Sinkhorn iteration). It is also possible to approximate transport using sums of Gaussians .

An interesting alternative is to compute a continuous function that approximates one of the pointsets (while keeping the other one discrete). This (assymetric) setting is thus called "semi-discrete", and can be solved numerically using efficient algorithms  (and my own one in ). It is implemented in my GEOGRAM software library .

 Optimal Transport Old and New - Cedric Villani

 Optimal Transport for Applied Mathematicians - Fillipo Santambrogio

 Displacement Interpolation Using Lagrangian Mass Transport, Bonneel et.al, ACM Transactions on Graphics (SIGGRAPH ASIA), 2011

 A multiscale approach to optimal transport, Quentin Merigot, Eurographics

 A numerical algorithm for semi-discrete optimal transport in 3D, Bruno Lévy, M2AN

• Thanks, Bruno. That looks like it might "work", i.e., accomplish a satisfactory-or-better visual appearance. Maybe even better-looking than my conjured-up procedure (or maybe not; I'd have to try both and take a look). As mentioned in comment to led23 above, I probably should have better-described my specific application, which would be like morphing a bitmap of the characters "Hello," into a bitmap of "World.", as an additional effect for my program forkosh.com/gifscroll.html Jun 16, 2017 at 22:18
• I'm finding the optimal transport problem itself very interesting, as per easy-to-follow cims.nyu.edu/~essid/Notes/Notes_optimal_transport.pdf and www-fourier.ujf-grenoble.fr/sites/ifmaquette.ujf-grenoble.fr/… in addition to your terrific references, including two free books. Code-wise, I found github.com/maxdan94/auction for the auction algorithm, and various stuff for the related hungarian algorithm, e.g., csclab.murraystate.edu/~bob.pilgrim/445/munkres.html . But I'm not finding much for sinkhorn. Do you know of good C code for that? Jun 17, 2017 at 8:14
• I think there is a Python implementation in Remy Flamary"s code: github.com/rflamary/POT (BTW, thank you for the additional references !). You may also lookup "earth movers distance" (that corresponds to $c(x,y) = \| x - y \|$). Jun 17, 2017 at 14:42
• Thanks again, Bruno. The .py stuff in POT-master/ot/ is a little (or >little) too dense for me to follow, especially given my next-to-nothing knowledge of python (and that's approaching zero from the negative side). Especially since I'd want to refactor it and translate it into C, in such a way that it fits neatly into the existing gifscroll.c code. All-in-all, led23's nearest-neighbor-interpolation seems the best compromise (perhaps heavily weighted towards "easy for me to do"). But the formalism of the optimal transport problem seems more and more interesting to me. Thanks so much for that. Jun 18, 2017 at 4:40