Optical flow works poorly on simulated data

Sample image and/or macro code

So I simulated simple data consisting of 2 matrices a and b. b is a deformed version of a. I want to find a transformation that will make b as a. Thus I decided on registration, and since I expect in my data not only translations, but also more nonlinear changes I decided on the optical flow algorithm. But the results of it look terrible. I don’t fully understand what OF does under the hood, but the scikit-image example from the docs (with the motorcycle) was quite impressive. Is there anything I’m doing wrong, or it wouldn’t work in this case? If so, why?

from skimage.transform import warp
from skimage.registration import optical_flow_tvl1
import numpy as np
import matplotlib.pyplot as plt

a = np.zeros((20, 20))
a[3:15,4:17] = 1
b = np.zeros((20, 20))
b[4:17,3:18] = 1
b[10:15,11:15] = 0

v, u = optical_flow_tvl1(a, b)

nr, nc = a.shape
r_coords, c_coords = np.meshgrid(np.arange(nr), np.arange(nc), indexing='ij')
b_warp = warp(mask1, np.array([r_coords + v, c_coords + u]), mode='nearest')
b_warp = (b_warp - np.min(b_warp)) / np.max(b_warp)

plt.figure(figsize=(15,3))
plt.subplot(131)
plt.imshow(a)
plt.subplot(132)
plt.imshow(b)
plt.subplot(133)
plt.imshow(b_warp)
plt.show()

From left to right you have a, b and b after OF transformation.

Analysis goals

Transform matrix b too look as close as possible to a.

Hi @dokato, welcome to the forum! This toy example is actually a very hard problem because the image has very little structure, its gradient is zero almost everywhere. If you want to simulate synthetic data for optical-flow registration you should use data with more structure. If you really need to register together these exact images please tell us a little more about what you would like to do :slight_smile:

Hello @dokato,

Quick answer, in your code you gave, the following line seems to be incorrect:
b_warp = warp(mask1, np.array([r_coords + v, c_coords + u]), mode='nearest')
when replaced with
b_warp = warp(b, np.array([r_coords + v, c_coords + u]), mode='nearest')
which I assume is what you want to do, it gives much better results:

That being said the result is not perfect but keep in mind that TV-L1, if I remember correctly, first has constraints on the continuity of the optic flow and second is supposed to work on images with continuous intensities and not binary images. These two constrains makes it hard to perfectly register your simulated data onto each other (as said @emmanuelle).

Now, if you make the resulting image binary simply with a glob threshold at .5 you get pretty good results:


(where the light blue is the mismatch between a and b_warp).

To generate the last figure I changed the following line:
plt.imshow(b_warp)
by that line:
plt.imshow((.5<b_warp)+a)

I hope this helps.

1 Like

The main hypothesis of TV-L1 optical flow are

  • the image intensity preservation: b[y+v, x+u] = a[y, x]
  • the differentiability of the inputs (a, b)
  • the differentiability of the solution (u, v).

@L.eo first catch a bug in your example and shows that the solution proposed by optical_flow_tvl1 is acceptable (b[y+v, x+u] ~ a[y, x]):

>>> a = np.zeros((20, 20))
>>> a[3:15,4:17] = 1
>>> b = np.zeros((20, 20))
>>> b[4:17,3:18] = 1
>>> b[10:15,11:15] = 0
>>> from skimage.registration import optical_flow_tvl1
>>> v, u = optical_flow_tvl1(a, b)
>>> from skimage.transform import warp
>>> r, c = np.meshgrid(np.arange(20), np.arange(20), indexing='ij')
>>> b_warp = warp(b, np.array([r+v, c+u]))

warp_bilinear
This can be highlighted by using order=0 in the warp function to apply nearest neighbor interpolation for image warping since a and b contains only two values:

>>> b_warp_2 = warp(b, np.array([r+v, c+u]), order=0)

warp_nn
The main problem of the toy problem that you implemented is that a large number of solution exist, and no one is able to say what is a correct solution!

Let us look at the shape of the solution proposed by optical_flow_tvl1:

>>> ax = plt.gca()
>>> ax = plt.gca()
>>> ax.axis([0, 20, 20, 0])
[0, 20, 20, 0]
>>> ax.set_aspect('equal')
>>> y, x = np.mgrid[:20, :20]
>>> P.quiver(x, y, u, v, angles='xy', units='dots', scale_units='xy')
<matplotlib.quiver.Quiver object at 0x7f40719f0690>
>>> P.imshow(a)
<matplotlib.image.AxesImage object at 0x7f4072357190>

quiver

Well, not that bad too ;)…