Example : Simple Lucas Kanade Optical Flow Estimator

In this example we present a simple approximation to the Lucas Kanade Optical Flow Estimator. The goal is to solve the 2x2 system describing the change in intensity by translation. We write the Pipeline Script LucasKanade.ppl which will contain both the algorithm implementation and a representation stage :

// The different formats needed :
TEXTURE_FORMAT:computationFormatFloat(512, 512,GL_RGBA32F, GL_FLOAT, GL_NEAREST, GL_NEAREST, GL_CLAMP, GL_CLAMP)
TEXTURE_FORMAT:visualizationFormat(512, 512, GL_RGB, GL_UNSIGNED_BYTE, GL_NEAREST, GL_NEAREST, GL_CLAMP, GL_CLAMP)
TEXTURE_FORMAT:sideBySideFormat(512, 256, GL_RGB, GL_UNSIGNED_BYTE, GL_NEAREST, GL_NEAREST, GL_CLAMP, GL_CLAMP)

SOURCE:derivativesShader
{
        #version 130

        uniform sampler2D       latest, oldest;
        out     vec4            derivatives, grayScale;

        void main()
        {
                const float step = 1.0/512.0;

                // Compute all the derivates of the intensity image :
                vec4 col11 = textureLod(latest, gl_TexCoord[0].st, 0);
                vec4 col01 = textureLod(latest, gl_TexCoord[0].st + vec2(-step, 0.0), 0);
                vec4 col21 = textureLod(latest, gl_TexCoord[0].st + vec2( step, 0.0), 0);
                vec4 col10 = textureLod(latest, gl_TexCoord[0].st + vec2(0.0, -step), 0);
                vec4 col12 = textureLod(latest, gl_TexCoord[0].st + vec2(0.0,  step), 0);
                vec4 colB  = textureLod(oldest, gl_TexCoord[0].st, 0);

                float   v11 = (col11.r+col11.g+col11.b)/3.0,
                        v01 = (col01.r+col01.g+col01.b)/3.0,
                        v21 = (col21.r+col21.g+col21.b)/3.0,
                        v10 = (col10.r+col10.g+col10.b)/3.0,
                        v12 = (col12.r+col12.g+col12.b)/3.0,
                        vB  = (colB.r+colB.g+colB.b)/3.0;

                derivatives.r = (v21-v01)/2.0;
                derivatives.g = (v12-v10)/2.0;
                derivatives.b = (v11-vB)/2.0;

                grayScale.rgb = vec3(v11,v11,v11);
        }
}

SOURCE:opticalFlowShader
{
        #version 130

        uniform sampler2D derivatives;
        out vec4 opticalFlow;

        void main()
        {
                const float step = 1.0/512.0;
                const int wSize = 5;
                const float mWindow = pow((2*wSize+1)*step/(3.14159265),2.0);
                float w = 0.0;
                int i, j;

                vec2    pos,
                        delta;
                vec4    col;
                float   x = 0.0,
                        y = 0.0,
                        A = 0.0,
                        B = 0.0,
                        C = 0.0,
                        idet = 0.0;

                // Build the 2x2 matrix and the Intensity vector: :
                for( i=-wSize; i<=wSize; i++)
                {
                        for( j=-wSize; j<=wSize; j++)
                        {
                                vec2 delta = vec2( i*step, j*step);
                                w       = cos(dot(delta,delta)*40.0);
                                pos     = gl_TexCoord[0].st + delta;
                                col     = textureLod(derivatives, pos, 0.0);
                                A       = A + w*col.r*col.r;
                                B       = B + w*col.g*col.g;
                                C       = C + w*col.r*col.g;
                                x       = x + w*col.r*col.b;
                                y       = y + w*col.g*col.b;
                        }
                }

                // Solve 2x2 matrix inverse :
                idet = 1.0/(A*B-C*C);

                if(idet>1000.0)
                        idet = 0.0;

                opticalFlow.r = idet*(-B*x+C*y);
                opticalFlow.g = idet*( C*x-A*y);
        }
}

SOURCE:visualizationShader
{
        #version 130

        uniform sampler2D opticalFlow;
        out vec4 visualization;

        const float eLim = 0.1;
        const float ctr = 0.1;

        void getHSL(in vec2 d, out vec3 res)
        {
                float   e       = (d.x*d.x + d.y*d.y)/2.0;

                // Angles :
                float   h       = (atan(d.x, d.y) + 3.14159)/1.047;
                int     prt     = int(h/2.0);
                float   hp      = h-float(prt)*2.0;
                float   v       = 1.0-abs(hp-1.0);

                if(h<1.0)
                        res = vec3(1.0,v,0.0)*e;
                else if(h<2.0)
                        res = vec3(v,1.0,0.0)*e;
                else if(h<3.0)
                        res = vec3(0.0,1.0,v)*e;
                else if(h<4.0)
                        res = vec3(0.0,v,1.0)*e;
                else if(h<5.0)
                        res = vec3(v,0.0,1.0)*e;
                else //if(h<6.0)
                        res = vec3(1.0,0.0,v)*e;
        }

        void main()
        {
                visualization.a = 1.0;

                // Draw a small legend for the directions :
                if(gl_TexCoord[0].s<ctr && gl_TexCoord[0].t<ctr)
                {
                        getHSL(gl_TexCoord[0].st-vec2(ctr,ctr)/2.0, visualization.rgb);
                        visualization.rgb = visualization.rgb*400.0;
                }
                else // The rest of the image :
                {
                        vec4 flow = textureLod(opticalFlow, gl_TexCoord[0].st, 0);

                        // Make the visualization easier by using HSL color space and associate directions with color.
                        getHSL(flow.rg, visualization.rgb);
                        visualization.rgb = visualization.rgb*0.5;
                }
        }
}

SOURCE:sideBySideShader
{
        #version 130

        uniform sampler2D latest, visualization;
        out vec4 sideBySide;

        void main()
        {
                vec2 pos = gl_TexCoord[0].st;

                // Separate the output image in left image and right image :
                if(pos.s<=0.5)
                {
                        // The original image goes on left :
                        pos.s = pos.s*2.0;
                        sideBySide  = textureLod(latest, pos, 0);
                }
                else
                {
                        // The motion image goes on right :
                        pos.s = (pos.s-0.5)*2.0;
                        sideBySide = textureLod(visualization, pos, 0);
                }
        }
}

FILTER_LAYOUT:derivatives(computationFormatFloat,derivativesShader)
FILTER_LAYOUT:opticalFlow(computationFormatFloat,opticalFlowShader)
FILTER_LAYOUT:visualizationFilter(visualizationFormat,visualizationShader)
FILTER_LAYOUT:sideBySideFilter(sideBySideFormat,sideBySideShader)

PIPELINE_MAIN:opticalFlow
{
        INPUT_PORTS(latest, oldest)
        OUTPUT_PORTS(grayScale, derivatives, opticalFlow, visualization, sideBySide)

        FILTER_INSTANCE:instDerivatives(cDerivatives)
        FILTER_INSTANCE:instOpticalFlow(cOpticalFlow)
        FILTER_INSTANCE:instVisualization(cMix)
        FILTER_INSTANCE:instSideBySide(cSBS)

        // Automatic connection will be performed.
}

Example of recording of the output port sideBySide with FFMPEG :

Animation Big Buck Bunny under Creative Commons License, see www.bigbuckbunny.org.