updated gpu optical_flow sample

modules/core/include/opencv2/core/opengl_interop.hpp

@@ -265,7 +265,7 @@ namespace cv
     }
 
     //! render OpenGL arrays
-    CV_EXPORTS void render(const GlArrays& arr, int mode = RenderMode::POINTS);
+    CV_EXPORTS void render(const GlArrays& arr, int mode = RenderMode::POINTS, Scalar color = Scalar::all(255));
 
     CV_EXPORTS void render(const std::string& str, const Ptr<GlFont>& font, Scalar color, Point2d pos);
 

modules/core/src/opengl_interop.cpp

@@ -1325,11 +1325,13 @@ void cv::render(const GlTexture& tex, Rect_<double> wndRect, Rect_<double> texRe
 #endif
 }
 
-void cv::render(const GlArrays& arr, int mode)
+void cv::render(const GlArrays& arr, int mode, Scalar color)
 {
 #ifndef HAVE_OPENGL
     throw_nogl;
 #else
+    glColor3d(color[0] / 255.0, color[1] / 255.0, color[3] / 255.0);
+
     arr.bind();
 
     glDrawArrays(mode, 0, arr.size().area());

modules/gpu/include/opencv2/gpu/gpu.hpp

@@ -1734,6 +1734,8 @@ CV_EXPORTS void interpolateFrames(const GpuMat& frame0, const GpuMat& frame1,
                                   float pos, GpuMat& newFrame, GpuMat& buf,
                                   Stream& stream = Stream::Null());
 
+CV_EXPORTS void createOpticalFlowNeedleMap(const GpuMat& u, const GpuMat& v, GpuMat& vertex, GpuMat& colors);
+
 } // namespace gpu
 
 //! Speckle filtering - filters small connected components on diparity image.

modules/gpu/src/cuda/canny.cu

@@ -81,7 +81,7 @@ namespace cv { namespace gpu { namespace device
             calcSobelRowPass<<<grid, block>>>(src, dx_buf, dy_buf, rows, cols);
             cudaSafeCall( cudaGetLastError() );
 
-            cudaSafeCall(cudaThreadSynchronize());
+            cudaSafeCall( cudaDeviceSynchronize() );
         }
 
         struct L1
@@ -171,7 +171,7 @@ namespace cv { namespace gpu { namespace device
 
             cudaSafeCall( cudaGetLastError() );
 
-            cudaSafeCall(cudaThreadSynchronize());
+            cudaSafeCall( cudaDeviceSynchronize() );
         }
 
         //////////////////////////////////////////////////////////////////////////////////////////
@@ -252,7 +252,7 @@ namespace cv { namespace gpu { namespace device
             calcMap<<<grid, block>>>(dx, dy, mag, map, rows, cols, low_thresh, high_thresh);
             cudaSafeCall( cudaGetLastError() );
 
-            cudaSafeCall(cudaThreadSynchronize());
+            cudaSafeCall( cudaDeviceSynchronize() );
         }
 
         //////////////////////////////////////////////////////////////////////////////////////////
@@ -345,7 +345,7 @@ namespace cv { namespace gpu { namespace device
             edgesHysteresisLocal<<<grid, block>>>(map, st1, rows, cols);
             cudaSafeCall( cudaGetLastError() );
 
-            cudaSafeCall(cudaThreadSynchronize());
+            cudaSafeCall( cudaDeviceSynchronize() );
         }
 
         __constant__ int c_dx[8] = {-1,  0,  1, -1, 1, -1, 0, 1};
@@ -460,7 +460,7 @@ namespace cv { namespace gpu { namespace device
                 edgesHysteresisGlobal<<<grid, block>>>(map, st1, st2, rows, cols, count);
                 cudaSafeCall( cudaGetLastError() );
 
-                cudaSafeCall(cudaThreadSynchronize());
+                cudaSafeCall( cudaDeviceSynchronize() );
 
                 cudaSafeCall( cudaMemcpy(&count, counter_ptr, sizeof(unsigned int), cudaMemcpyDeviceToHost) );
 
@@ -485,7 +485,7 @@ namespace cv { namespace gpu { namespace device
             getEdges<<<grid, block>>>(map, dst, rows, cols);
             cudaSafeCall( cudaGetLastError() );
 
-            cudaSafeCall(cudaThreadSynchronize());
+            cudaSafeCall( cudaDeviceSynchronize() );
         }
     } // namespace canny
 }}} // namespace cv { namespace gpu { namespace device

modules/gpu/src/cuda/optical_flow.cu

@@ -0,0 +1,217 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "opencv2/gpu/device/common.hpp"
+
+namespace cv { namespace gpu { namespace device 
+{
+    namespace optical_flow
+    {
+        #define NEEDLE_MAP_SCALE 16
+        #define MAX_FLOW 30.0f
+        #define NUM_VERTS_PER_ARROW 6
+
+        __global__ void NeedleMapAverageKernel(const DevMem2Df u, const PtrStepf v, PtrStepf u_avg, PtrStepf v_avg)
+        {   
+            __shared__ float smem[2 * NEEDLE_MAP_SCALE];
+
+            volatile float* u_col_sum = smem;
+            volatile float* v_col_sum = u_col_sum + NEEDLE_MAP_SCALE;
+
+            const int x = blockIdx.x * NEEDLE_MAP_SCALE + threadIdx.x;
+            const int y = blockIdx.y * NEEDLE_MAP_SCALE;
+
+            u_col_sum[threadIdx.x] = 0;
+            v_col_sum[threadIdx.x] = 0;
+
+            #pragma unroll
+            for(int i = 0; i < NEEDLE_MAP_SCALE; ++i)
+            {
+                u_col_sum[threadIdx.x] += u(::min(y + i, u.rows - 1), x);
+                v_col_sum[threadIdx.x] += v(::min(y + i, u.rows - 1), x);
+            }
+
+            if (threadIdx.x < 8)
+            {        	
+                // now add the column sums
+                const uint X = threadIdx.x;
+
+                if (X | 0xfe == 0xfe)  // bit 0 is 0
+                {
+                    u_col_sum[threadIdx.x] += u_col_sum[threadIdx.x + 1];
+                    v_col_sum[threadIdx.x] += v_col_sum[threadIdx.x + 1];
+                }
+
+                if (X | 0xfe == 0xfc) // bits 0 & 1 == 0 
+                { 
+                    u_col_sum[threadIdx.x] += u_col_sum[threadIdx.x + 2];
+                    v_col_sum[threadIdx.x] += v_col_sum[threadIdx.x + 2];
+                }
+
+                if (X | 0xf8 == 0xf8)
+                {
+                    u_col_sum[threadIdx.x] += u_col_sum[threadIdx.x + 4];
+                    v_col_sum[threadIdx.x] += v_col_sum[threadIdx.x + 4];
+                }
+
+                if (X == 0)
+                {
+                    u_col_sum[threadIdx.x] += u_col_sum[threadIdx.x + 8];
+                    v_col_sum[threadIdx.x] += v_col_sum[threadIdx.x + 8];
+                }
+            }
+
+            if (threadIdx.x == 0)
+            {
+                const float coeff = 1.0f / (NEEDLE_MAP_SCALE * NEEDLE_MAP_SCALE);
+
+                u_col_sum[0] *= coeff;
+                v_col_sum[0] *= coeff;
+
+                u_avg(blockIdx.y, blockIdx.x) = u_col_sum[0];
+                v_avg(blockIdx.y, blockIdx.x) = v_col_sum[0];
+            }
+        }
+        
+        void NeedleMapAverage_gpu(DevMem2Df u, DevMem2Df v, DevMem2Df u_avg, DevMem2Df v_avg)
+        {
+            const dim3 block(NEEDLE_MAP_SCALE);
+            const dim3 grid(u_avg.cols, u_avg.rows);
+
+            NeedleMapAverageKernel<<<grid, block>>>(u, v, u_avg, v_avg);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+
+        __global__ void NeedleMapVertexKernel(const DevMem2Df u_avg, const PtrStepf v_avg, float* vertex_data, float* color_data, float xscale, float yscale)
+        {
+            // test - just draw a triangle at each pixel
+            const int x = blockIdx.x * blockDim.x + threadIdx.x;
+            const int y = blockIdx.y * blockDim.y + threadIdx.y;
+
+            const float arrow_x = x * NEEDLE_MAP_SCALE + NEEDLE_MAP_SCALE / 2.0f;
+            const float arrow_y = y * NEEDLE_MAP_SCALE + NEEDLE_MAP_SCALE / 2.0f;
+
+            float3 v[NUM_VERTS_PER_ARROW];
+
+            if (x < u_avg.cols && y < u_avg.rows)
+            {
+                const float u_avg_val = u_avg(y, x);
+                const float v_avg_val = v_avg(y, x);
+
+                const float theta = ::atan2f(v_avg_val, u_avg_val) + CV_PI;
+
+                float r = ::sqrtf(v_avg_val * v_avg_val + u_avg_val * u_avg_val);
+                r = fmin(14.0f * (r / MAX_FLOW), 14.0f);
+
+                v[0].z = 1.0f;
+                v[1].z = 0.7f;
+                v[2].z = 0.7f;
+                v[3].z = 0.7f;
+                v[4].z = 0.7f;
+                v[5].z = 1.0f;
+
+                v[0].x = arrow_x;
+                v[0].y = arrow_y;
+                v[5].x = arrow_x;
+                v[5].y = arrow_y;
+
+                v[2].x = arrow_x + r * ::cosf(theta);
+                v[2].y = arrow_y + r * ::sinf(theta);
+                v[3].x = v[2].x;
+                v[3].y = v[2].y;
+
+                r = ::fmin(r, 2.5f);
+
+                v[1].x = arrow_x + r * ::cosf(theta - CV_PI / 2.0f);
+                v[1].y = arrow_y + r * ::sinf(theta - CV_PI / 2.0f);
+
+                v[4].x = arrow_x + r * ::cosf(theta + CV_PI / 2.0f);
+                v[4].y = arrow_y + r * ::sinf(theta + CV_PI / 2.0f);
+
+                int indx = (y * u_avg.cols + x) * NUM_VERTS_PER_ARROW * 3;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[0].x * xscale;
+                vertex_data[indx++] = v[0].y * yscale;
+                vertex_data[indx++] = v[0].z;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[1].x * xscale;
+                vertex_data[indx++] = v[1].y * yscale;
+                vertex_data[indx++] = v[1].z;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[2].x * xscale;
+                vertex_data[indx++] = v[2].y * yscale;
+                vertex_data[indx++] = v[2].z;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[3].x * xscale;
+                vertex_data[indx++] = v[3].y * yscale;
+                vertex_data[indx++] = v[3].z;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[4].x * xscale;
+                vertex_data[indx++] = v[4].y * yscale;
+                vertex_data[indx++] = v[4].z;
+
+                color_data[indx] = (theta - CV_PI) / CV_PI * 180.0f;
+                vertex_data[indx++] = v[5].x * xscale;
+                vertex_data[indx++] = v[5].y * yscale;
+                vertex_data[indx++] = v[5].z;
+            }
+        }
+
+        void CreateOpticalFlowNeedleMap_gpu(DevMem2Df u_avg, DevMem2Df v_avg, float* vertex_buffer, float* color_data, float xscale, float yscale)
+        {
+            const dim3 block(16);
+            const dim3 grid(divUp(u_avg.cols, block.x), divUp(u_avg.rows, block.y));
+
+            NeedleMapVertexKernel<<<grid, block>>>(u_avg, v_avg, vertex_buffer, color_data, xscale, yscale);
+            cudaSafeCall( cudaGetLastError() );
+
+            cudaSafeCall( cudaDeviceSynchronize() );
+        }
+    }
+}}}

modules/gpu/src/optical_flow.cpp

@@ -50,6 +50,7 @@ using namespace std;
 
 void cv::gpu::BroxOpticalFlow::operator ()(const GpuMat&, const GpuMat&, GpuMat&, GpuMat&, Stream&) { throw_nogpu(); }
 void cv::gpu::interpolateFrames(const GpuMat&, const GpuMat&, const GpuMat&, const GpuMat&, const GpuMat&, const GpuMat&, float, GpuMat&, GpuMat&, Stream&) { throw_nogpu(); }
+void cv::gpu::createOpticalFlowNeedleMap(const GpuMat&, const GpuMat&, GpuMat&, GpuMat&) { throw_nogpu(); }
 
 #else
 
@@ -188,4 +189,44 @@ void cv::gpu::interpolateFrames(const GpuMat& frame0, const GpuMat& frame1, cons
         cudaSafeCall( cudaDeviceSynchronize() );
 }
 
+namespace cv { namespace gpu { namespace device 
+{
+    namespace optical_flow
+    {
+        void NeedleMapAverage_gpu(DevMem2Df u, DevMem2Df v, DevMem2Df u_avg, DevMem2Df v_avg);
+        void CreateOpticalFlowNeedleMap_gpu(DevMem2Df u_avg, DevMem2Df v_avg, float* vertex_buffer, float* color_data, float xscale, float yscale);
+    }
+}}}
+
+void cv::gpu::createOpticalFlowNeedleMap(const GpuMat& u, const GpuMat& v, GpuMat& vertex, GpuMat& colors)
+{
+    using namespace cv::gpu::device::optical_flow;
+
+    CV_Assert(u.type() == CV_32FC1);
+    CV_Assert(v.type() == u.type() && v.size() == u.size());
+
+    const int NEEDLE_MAP_SCALE = 16;
+
+	const int x_needles = u.cols / NEEDLE_MAP_SCALE;
+	const int y_needles = u.rows / NEEDLE_MAP_SCALE;
+
+    GpuMat u_avg(y_needles, x_needles, CV_32FC1);
+    GpuMat v_avg(y_needles, x_needles, CV_32FC1);
+    
+    NeedleMapAverage_gpu(u, v, u_avg, v_avg);
+    
+    const int NUM_VERTS_PER_ARROW = 6;
+    
+    const int num_arrows = x_needles * y_needles * NUM_VERTS_PER_ARROW;
+
+    vertex.create(1, num_arrows, CV_32FC3);
+    colors.create(1, num_arrows, CV_32FC3);
+
+    colors.setTo(Scalar::all(1.0));
+
+    CreateOpticalFlowNeedleMap_gpu(u_avg, v_avg, vertex.ptr<float>(), colors.ptr<float>(), 1.0f / u.cols, 1.0f / u.rows);
+
+    cvtColor(colors, colors, COLOR_HSV2RGB);
+}
+
 #endif /* HAVE_CUDA */

samples/gpu/optical_flow.cpp

@@ -4,132 +4,245 @@
 
 #include "cvconfig.h"
 #include "opencv2/core/core.hpp"
+#include "opencv2/core/opengl_interop.hpp"
 #include "opencv2/highgui/highgui.hpp"
 #include "opencv2/gpu/gpu.hpp"
 
-#ifdef HAVE_CUDA
-#include "NPP_staging/NPP_staging.hpp"
-#endif
-
 using namespace std;
 using namespace cv;
 using namespace cv::gpu;
 
-#if !defined(HAVE_CUDA)
+void getFlowField(const Mat& u, const Mat& v, Mat& flowField);
+
+#ifdef HAVE_OPENGL
+
+void needleMapDraw(void* userdata);
+
+#endif
 
 int main(int argc, const char* argv[])
 {
-    cout << "Please compile the library with CUDA support" << endl;
+    try
-    return -1;
+    {
-}
+        const char* keys =
+           "{ h  | help      | false | print help message }"
+           "{ l  | left      |       | specify left image }"
+           "{ r  | right     |       | specify right image }"
+           "{ s  | scale     | 0.8   | set pyramid scale factor }"
+           "{ a  | alpha     | 0.197 | set alpha }"
+           "{ g  | gamma     | 50.0  | set gamma }"
+           "{ i  | inner     | 10    | set number of inner iterations }"
+           "{ o  | outer     | 77    | set number of outer iterations }"
+           "{ si | solver    | 10    | set number of basic solver iterations }"
+           "{ t  | time_step | 0.1   | set frame interpolation time step }";
 
-#else
+        CommandLineParser cmd(argc, argv, keys);
 
-#define PARAM_LEFT      "--left"
+        if (cmd.get<bool>("help"))
-#define PARAM_RIGHT     "--right"
-#define PARAM_SCALE     "--scale"
-#define PARAM_ALPHA     "--alpha"
-#define PARAM_GAMMA     "--gamma"
-#define PARAM_INNER     "--inner"
-#define PARAM_OUTER     "--outer"
-#define PARAM_SOLVER    "--solver"
-#define PARAM_TIME_STEP "--time_step"
-#define PARAM_HELP      "--help"
-
-bool help_showed = false;
-
-void printHelp()
         {
-    cout << "Usage help:\n";
+            cout << "Usage: optical_float [options]" << endl;
-    cout << setiosflags(ios::left);
+            cout << "Avaible options:" << endl;
-    cout << "\t" << setw(15) << PARAM_ALPHA << " - set alpha\n";
+            cmd.printParams();
-    cout << "\t" << setw(15) << PARAM_GAMMA << " - set gamma\n";
-    cout << "\t" << setw(15) << PARAM_INNER << " - set number of inner iterations\n";
-    cout << "\t" << setw(15) << PARAM_LEFT << " - specify left image\n";
-    cout << "\t" << setw(15) << PARAM_RIGHT << " - specify right image\n";
-    cout << "\t" << setw(15) << PARAM_OUTER << " - set number of outer iterations\n";
-    cout << "\t" << setw(15) << PARAM_SCALE << " - set pyramid scale factor\n";
-    cout << "\t" << setw(15) << PARAM_SOLVER << " - set number of basic solver iterations\n";
-    cout << "\t" << setw(15) << PARAM_TIME_STEP << " - set frame interpolation time step\n";
-    cout << "\t" << setw(15) << PARAM_HELP << " - display this help message\n";
-    help_showed = true;
-}
-
-int processCommandLine(int argc, const char* argv[], float& timeStep, string& frame0Name, string& frame1Name, BroxOpticalFlow& flow)
-{
-    timeStep = 0.25f;
-
-    for (int iarg = 1; iarg < argc; ++iarg)
-    {
-        if (strcmp(argv[iarg], PARAM_LEFT) == 0)
-        {
-            if (iarg + 1 < argc)
-                frame0Name = argv[++iarg];
-            else
-                return -1;
-        }
-        if (strcmp(argv[iarg], PARAM_RIGHT) == 0)
-        {
-            if (iarg + 1 < argc)
-                frame1Name = argv[++iarg];
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_SCALE) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.scale_factor = static_cast<float>(atof(argv[++iarg]));
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_ALPHA) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.alpha = static_cast<float>(atof(argv[++iarg]));
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_GAMMA) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.gamma = static_cast<float>(atof(argv[++iarg]));
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_INNER) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.inner_iterations = atoi(argv[++iarg]);
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_OUTER) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.outer_iterations = atoi(argv[++iarg]);
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_SOLVER) == 0)
-        {
-            if (iarg + 1 < argc)
-                flow.solver_iterations = atoi(argv[++iarg]);
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_TIME_STEP) == 0)
-        {
-            if (iarg + 1 < argc)
-                timeStep = static_cast<float>(atof(argv[++iarg]));
-            else
-                return -1;
-        }
-        else if(strcmp(argv[iarg], PARAM_HELP) == 0)
-        {
-            printHelp();
             return 0;
         }
+
+        string frame0Name = cmd.get<string>("left");
+        string frame1Name = cmd.get<string>("right");
+        float scale = cmd.get<float>("scale");
+        float alpha = cmd.get<float>("alpha");
+        float gamma = cmd.get<float>("gamma");
+        int inner_iterations = cmd.get<int>("inner");
+        int outer_iterations = cmd.get<int>("outer");
+        int solver_iterations = cmd.get<int>("solver");
+        float timeStep = cmd.get<float>("time_step");
+
+        if (frame0Name.empty() || frame1Name.empty())
+        {
+            cerr << "Missing input file names" << endl;
+            return -1;
         }
+
+        Mat frame0Color = imread(frame0Name);
+        Mat frame1Color = imread(frame1Name);
+
+        if (frame0Color.empty() || frame1Color.empty())
+        {
+            cout << "Can't load input images" << endl;
+            return -1;
+        }
+
+        cout << "OpenCV / NVIDIA Computer Vision" << endl;
+        cout << "Optical Flow Demo: Frame Interpolation" << endl;
+        cout << "=========================================" << endl;
+
+        namedWindow("Forward flow");
+        namedWindow("Backward flow");
+
+        namedWindow("Needle Map", WINDOW_OPENGL);
+
+        namedWindow("Interpolated frame");
+
+        setGlDevice();
+
+        cout << "Press:" << endl;
+        cout << "\tESC to quit" << endl;
+        cout << "\t'a' to move to the previous frame" << endl;
+        cout << "\t's' to move to the next frame\n" << endl;
+
+        frame0Color.convertTo(frame0Color, CV_32F, 1.0 / 255.0);
+        frame1Color.convertTo(frame1Color, CV_32F, 1.0 / 255.0);
+
+        Mat frame0Gray, frame1Gray;
+
+        cvtColor(frame0Color, frame0Gray, COLOR_BGR2GRAY);
+        cvtColor(frame1Color, frame1Gray, COLOR_BGR2GRAY);
+
+        GpuMat d_frame0(frame0Gray);
+        GpuMat d_frame1(frame1Gray);
+
+        cout << "Estimating optical flow" << endl;
+
+        BroxOpticalFlow d_flow(alpha, gamma, scale, inner_iterations, outer_iterations, solver_iterations);
+
+        cout << "\tForward..." << endl;
+
+        GpuMat d_fu, d_fv;
+
+        d_flow(d_frame0, d_frame1, d_fu, d_fv);
+                
+        Mat flowFieldForward;
+        getFlowField(Mat(d_fu), Mat(d_fv), flowFieldForward);
+    
+        cout << "\tBackward..." << endl;
+
+        GpuMat d_bu, d_bv;
+
+        d_flow(d_frame1, d_frame0, d_bu, d_bv);
+        
+        Mat flowFieldBackward;
+        getFlowField(Mat(d_bu), Mat(d_bv), flowFieldBackward);
+
+#ifdef HAVE_OPENGL
+        cout << "Create Optical Flow Needle Map..." << endl;
+        
+        GpuMat d_vertex, d_colors;
+
+        createOpticalFlowNeedleMap(d_bu, d_bv, d_vertex, d_colors);
+#endif
+
+        cout << "Interpolating..." << endl;
+
+        // first frame color components
+        GpuMat d_b, d_g, d_r;
+
+        // second frame color components
+        GpuMat d_bt, d_gt, d_rt;
+
+        // prepare color components on host and copy them to device memory
+        Mat channels[3];
+        cv::split(frame0Color, channels);
+
+        d_b.upload(channels[0]);
+        d_g.upload(channels[1]);
+        d_r.upload(channels[2]);
+
+        cv::split(frame1Color, channels);
+
+        d_bt.upload(channels[0]);
+        d_gt.upload(channels[1]);
+        d_rt.upload(channels[2]);
+    
+        // temporary buffer
+        GpuMat d_buf;
+
+        // intermediate frame color components (GPU memory)
+        GpuMat d_rNew, d_gNew, d_bNew;
+
+        GpuMat d_newFrame;
+
+        vector<Mat> frames;
+        frames.reserve(static_cast<int>(1.0f / timeStep) + 2);
+
+        frames.push_back(frame0Color);
+
+        // compute interpolated frames
+        for (float timePos = timeStep; timePos < 1.0f; timePos += timeStep)
+        {
+            // interpolate blue channel
+            interpolateFrames(d_b, d_bt, d_fu, d_fv, d_bu, d_bv, timePos, d_bNew, d_buf);
+
+            // interpolate green channel
+            interpolateFrames(d_g, d_gt, d_fu, d_fv, d_bu, d_bv, timePos, d_gNew, d_buf);
+
+            // interpolate red channel
+            interpolateFrames(d_r, d_rt, d_fu, d_fv, d_bu, d_bv, timePos, d_rNew, d_buf);
+
+            GpuMat channels[] = {d_bNew, d_gNew, d_rNew};
+            merge(channels, 3, d_newFrame);
+
+            frames.push_back(Mat(d_newFrame));
+
+            cout << setprecision(4) << timePos * 100.0f << "%\r";
+        }
+
+        frames.push_back(frame1Color);
+
+        cout << setw(5) << "100%" << endl;
+
+        cout << "Done" << endl;
+
+        imshow("Forward flow", flowFieldForward);
+        imshow("Backward flow", flowFieldBackward);
+
+#ifdef HAVE_OPENGL
+        GlArrays arr;
+        arr.setVertexArray(d_vertex);
+        arr.setColorArray(d_colors, false);
+
+        setOpenGlDrawCallback("Needle Map", needleMapDraw, &arr);
+#endif
+
+        int currentFrame = 0;
+
+        imshow("Interpolated frame", frames[currentFrame]);
+
+        while (true)
+        {
+            int key = toupper(waitKey(10));
+
+            switch (key)
+            {
+            case 27:
+                return 0;
+                break;
+
+            case 'A':
+                if (currentFrame > 0) 
+                    --currentFrame;
+
+                imshow("Interpolated frame", frames[currentFrame]);
+                break;
+
+            case 'S':
+                if (currentFrame < frames.size() - 1)
+                    ++currentFrame;
+
+                imshow("Interpolated frame", frames[currentFrame]);
+                break;
+            }
+        }
+    }
+    catch (const exception& ex)
+    {
+        cerr << ex.what() << endl;
+        return -1;
+    }
+    catch (...)
+    {
+        cerr << "Unknow error" << endl;
+        return -1;
+    }
+
     return 0;
 }
 
@@ -181,171 +294,20 @@ void getFlowField(const Mat& u, const Mat& v, Mat& flowField)
     }
 }
 
-int main(int argc, const char* argv[])
+#ifdef HAVE_OPENGL
+
+void needleMapDraw(void* userdata)
 {
-    string frame0Name, frame1Name;
+    const GlArrays* arr = static_cast<const GlArrays*>(userdata);
-    float timeStep = 0.01f;
 
-    BroxOpticalFlow d_flow(0.197f /*alpha*/, 50.0f /*gamma*/, 0.8f /*scale_factor*/, 
+    GlCamera camera;
-                           10 /*inner_iterations*/, 77 /*outer_iterations*/, 10 /*solver_iterations*/);
+    camera.setOrthoProjection(0.0, 1.0, 1.0, 0.0, 0.0, 1.0);
+    camera.lookAt(Point3d(0.0, 0.0, 1.0), Point3d(0.0, 0.0, 0.0), Point3d(0.0, 1.0, 0.0));
 
-    int result = processCommandLine(argc, argv, timeStep, frame0Name, frame1Name, d_flow);
+    camera.setupProjectionMatrix();
-    if (help_showed)
+    camera.setupModelViewMatrix();
-        return -1;
-    if (argc == 1 || result)
-    {
-        printHelp();
-        return result;
-    }
 
-    if (frame0Name.empty() || frame1Name.empty())
+    render(*arr, RenderMode::TRIANGLES);
-    {
-        cout << "Missing input file names\n";
-        return -1;
-    }
-
-    Mat frame0Color = imread(frame0Name);
-    Mat frame1Color = imread(frame1Name);
-
-    if (frame0Color.empty() || frame1Color.empty())
-    {
-        cout << "Can't load input images\n";
-        return -1;
-    }
-
-    cout << "OpenCV / NVIDIA Computer Vision\n";
-    cout << "Optical Flow Demo: Frame Interpolation\n";
-    cout << "=========================================\n";
-    cout << "Press:\n ESC to quit\n 'a' to move to the previous frame\n 's' to move to the next frame\n";
-
-    frame0Color.convertTo(frame0Color, CV_32F, 1.0 / 255.0);
-    frame1Color.convertTo(frame1Color, CV_32F, 1.0 / 255.0);
-
-    Mat frame0Gray, frame1Gray;
-
-    cvtColor(frame0Color, frame0Gray, COLOR_BGR2GRAY);
-    cvtColor(frame1Color, frame1Gray, COLOR_BGR2GRAY);
-
-    GpuMat d_frame0(frame0Gray);
-    GpuMat d_frame1(frame1Gray);
-    
-    Mat fu, fv;
-    Mat bu, bv;
-
-    GpuMat d_fu, d_fv;
-    GpuMat d_bu, d_bv;
-
-    cout << "Estimating optical flow\nForward...\n";
-
-    d_flow(d_frame0, d_frame1, d_fu, d_fv);
-    d_flow(d_frame1, d_frame0, d_bu, d_bv);
-    
-    d_fu.download(fu);
-    d_fv.download(fv);
-    
-    d_bu.download(bu);
-    d_bv.download(bv);
-
-    // first frame color components (GPU memory)
-    GpuMat d_b, d_g, d_r;
-
-    // second frame color components (GPU memory)
-    GpuMat d_bt, d_gt, d_rt;
-
-    // prepare color components on host and copy them to device memory
-    Mat channels[3];
-
-    cv::split(frame0Color, channels);
-
-    d_b.upload(channels[0]);
-    d_g.upload(channels[1]);
-    d_r.upload(channels[2]);
-
-    cv::split(frame1Color, channels);
-
-    d_bt.upload(channels[0]);
-    d_gt.upload(channels[1]);
-    d_rt.upload(channels[2]);
-
-    cout << "Interpolating...\n";
-    cout.precision (4);
-    
-    // temporary buffer
-    GpuMat d_buf;
-
-    // intermediate frame color components (GPU memory)
-    GpuMat d_rNew, d_gNew, d_bNew;
-
-    GpuMat d_newFrame;
-
-    vector<Mat> frames;
-    frames.reserve(1.0f / timeStep + 2);
-
-    frames.push_back(frame0Color);
-
-    // compute interpolated frames
-    for (float timePos = timeStep; timePos < 1.0f; timePos += timeStep)
-    {
-        // interpolate blue channel
-        interpolateFrames(d_b, d_bt, d_fu, d_fv, d_bu, d_bv, timePos, d_bNew, d_buf);
-        // interpolate green channel
-        interpolateFrames(d_g, d_gt, d_fu, d_fv, d_bu, d_bv, timePos, d_gNew, d_buf);
-        // interpolate red channel
-        interpolateFrames(d_r, d_rt, d_fu, d_fv, d_bu, d_bv, timePos, d_rNew, d_buf);
-
-        GpuMat channels[] = {d_bNew, d_gNew, d_rNew};
-        merge(channels, 3, d_newFrame);
-
-        Mat newFrame;
-        d_newFrame.download(newFrame);
-
-        frames.push_back(newFrame);
-
-        cout << timePos * 100.0f << "%\r";
-    }
-    cout << setw (5) << "100%\n";
-
-    frames.push_back(frame1Color);
-
-    int currentFrame;
-    currentFrame = 0;
-
-    Mat flowFieldForward;
-    Mat flowFieldBackward;
-
-    getFlowField(fu, fv, flowFieldForward);
-    getFlowField(bu, bv, flowFieldBackward);
-
-    imshow("Forward flow", flowFieldForward);
-    imshow("Backward flow", flowFieldBackward);
-
-    imshow("Interpolated frame", frames[currentFrame]);
-
-    bool qPressed = false;
-    while (!qPressed)
-    {
-        int key = toupper(waitKey(10));
-        switch (key)
-        {
-        case 27:
-            qPressed = true;
-            break;
-        case 'A':
-            if (currentFrame > 0) 
-                --currentFrame;
-
-            imshow("Interpolated frame", frames[currentFrame]);
-            break;
-        case 'S':
-            if (currentFrame < frames.size() - 1)
-                ++currentFrame;
-
-            imshow("Interpolated frame", frames[currentFrame]);
-            break;
-        }
-    }
-
-    return 0;
 }
 
 #endif