# 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

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]))
``````

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)
``````

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>
``````

Well, not that bad too ;)…