Implementing Flow operator

[carterbox]

This PR is supposed to introduce the Flow operator and wrap Farnebacks flow solver from OpenCV, but I am having trouble getting Farnebacks method to work for simple uniform shifts!

UPDATE: I discovered that the OpenCV implementation is sensitive to receiving view of arrays. So instead of accepting any shaped inputs, the Farneback alignment method can only do 3D arrays. I created some tests for the flow operator and Farneback method. The Flow operator's inverse is only adjoint for integer shifts. This may be due to interpolation errors.

[pep8speaks]

Hello @carterbox! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! :beers:

Comment last updated at 2020-07-08 22:54:53 UTC

[nikitinvv]

@carterbox, input arrays for Farneback's algorithm have different amplitudes whenever you convert them to uint8,

You should use the same range for normalization,

a = np.angle(data[0])
b = np.angle(original[0])
a=np.uint8((a-np.min(b))/(np.max(b)-np.min(b))*255)
b=np.uint8((b-np.min(b))/(np.max(b)-np.min(b))*255)

s = cv2.calcOpticalFlowFarneback(
            a,b,
            flow=None,
            pyr_scale=0.5,
            levels=3,
            winsize=12,
            iterations=4,
            poly_n=5,
            poly_sigma=1.1, 
            flags=0,
)
print(a.shape)
plt.subplot(1,2,1)
plt.imshow(a)
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(b)
plt.colorbar()
print(s[64,64,:])

Now, I would determine the range by using a histogram for unaligned data at the beginning. The histogram allows dealing with noisy data. For example, if in your image you put one pixel as a very huge value, then your approach with normalization to [0,255] with computing min,max for 2 images simply fails.

This is how I would determine the range [mmin,mmax]

def find_min_max(data):
    """Find min and max values according to histogram"""

    mmin = np.zeros(data.shape[0],dtype='float32')
    mmax = np.zeros(data.shape[0],dtype='float32')

    for k in range(data.shape[0]):
        h, e = np.histogram(data[k][:],1000)
        stend = np.where(h>np.max(h)*0.005)
        st = stend[0][0]
        end = stend[0][-1]        
        mmin[k] = e[st]
        mmax[k] = e[end+1]

    return mmin,mmax

Then I would perform normalization to [0,255] according to mmin,mmax with cutting values that are out of this interval. Something like this

  def registration_flow(self, psi, g, mmin, mmax, flow, pars, id):
        """Find optical flow for one projection"""
        tmp1 = ((psi[id]-mmin[id]) /
                (mmax[id]-mmin[id])*255)
        tmp1[tmp1 > 255] = 255
        tmp1[tmp1 < 0] = 0
        tmp2 = ((g[id]-mmin[id]) /
                (mmax[id]-mmin[id])*255)
        tmp2[tmp2 > 255] = 255
        tmp2[tmp2 < 0] = 0
        cv2.calcOpticalFlowFarneback(
           tmp1, tmp2, flow[id], *pars)  # updates flow

good luck

[carterbox]

I have determined that there is a problem with OpenCV. If the input arrays are reshaped before being passed to calcOpticalFlowFarneback, even if they don't actually change shape, OpenCV returns garbage.

[nikitinvv]

@carterbox np.ascontiguousarray doesnt help?

[carterbox]

@nikitinvv No, ascontiguousarray doesn't help.

[carterbox]

@nikitinvv I guess that would explain why I was having trouble with initial guesses! Do you know where I can find documentation for those flags? I couldn't find any.

[carterbox]

Could you add some check for this?

What kind of check were you thinking?

[nikitinvv]

@carterbox, I spent quite some time to find this 4 :)). Finally, I looked directly into C++ farnebacks code in opencv and found some #define ... 4 statements in a header. It is useful to have a code like ' if flow is not none and flag==4' then error or give a warning about that the flag will not work. Maybe lets even not allow users to choose this flag since it is kindly strange. We could do: If flow is none then flag = 0, else flag=4.