modules/cudaoptflow/src/tvl1flow.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44
  45 #if !defined HAVE_CUDA || defined(CUDA_DISABLER)
  46
  47 cv::cuda::OpticalFlowDual_TVL1_CUDA::OpticalFlowDual_TVL1_CUDA() { throw_no_cuda(); }
  48 void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat&, const GpuMat&, GpuMat&, GpuMat&) { throw_no_cuda(); }
  49 void cv::cuda::OpticalFlowDual_TVL1_CUDA::collectGarbage() {}
  50 void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat&, const GpuMat&, GpuMat&, GpuMat&) { throw_no_cuda(); }
  51
  52 #else
  53
  54 using namespace cv;
  55 using namespace cv::cuda;
  56
  57 cv::cuda::OpticalFlowDual_TVL1_CUDA::OpticalFlowDual_TVL1_CUDA()
  58 {
  59     tau            = 0.25;
  60     lambda         = 0.15;
  61     theta          = 0.3;
  62     nscales        = 5;
  63     warps          = 5;
  64     epsilon        = 0.01;
  65     iterations     = 300;
  66     scaleStep      = 0.8;
  67     gamma           = 0.0;
  68     useInitialFlow = false;
  69 }
  70
  71 void cv::cuda::OpticalFlowDual_TVL1_CUDA::operator ()(const GpuMat& I0, const GpuMat& I1, GpuMat& flowx, GpuMat& flowy)
  72 {
  73     CV_Assert( I0.type() == CV_8UC1 || I0.type() == CV_32FC1 );
  74     CV_Assert( I0.size() == I1.size() );
  75     CV_Assert( I0.type() == I1.type() );
  76     CV_Assert( !useInitialFlow || (flowx.size() == I0.size() && flowx.type() == CV_32FC1 && flowy.size() == flowx.size() && flowy.type() == flowx.type()) );
  77     CV_Assert( nscales > 0 );
  78
  79     // allocate memory for the pyramid structure
  80     I0s.resize(nscales);
  81     I1s.resize(nscales);
  82     u1s.resize(nscales);
  83     u2s.resize(nscales);
  84     u3s.resize(nscales);
  85
  86     I0.convertTo(I0s[0], CV_32F, I0.depth() == CV_8U ? 1.0 : 255.0);
  87     I1.convertTo(I1s[0], CV_32F, I1.depth() == CV_8U ? 1.0 : 255.0);
  88
  89     if (!useInitialFlow)
  90     {
  91         flowx.create(I0.size(), CV_32FC1);
  92         flowy.create(I0.size(), CV_32FC1);
  93     }
  94
  95     u1s[0] = flowx;
  96     u2s[0] = flowy;
  97     if (gamma)
  98         u3s[0].create(I0.size(), CV_32FC1);
  99
 100     I1x_buf.create(I0.size(), CV_32FC1);
 101     I1y_buf.create(I0.size(), CV_32FC1);
 102
 103     I1w_buf.create(I0.size(), CV_32FC1);
 104     I1wx_buf.create(I0.size(), CV_32FC1);
 105     I1wy_buf.create(I0.size(), CV_32FC1);
 106
 107     grad_buf.create(I0.size(), CV_32FC1);
 108     rho_c_buf.create(I0.size(), CV_32FC1);
 109
 110     p11_buf.create(I0.size(), CV_32FC1);
 111     p12_buf.create(I0.size(), CV_32FC1);
 112     p21_buf.create(I0.size(), CV_32FC1);
 113     p22_buf.create(I0.size(), CV_32FC1);
 114     if (gamma)
 115     {
 116         p31_buf.create(I0.size(), CV_32FC1);
 117         p32_buf.create(I0.size(), CV_32FC1);
 118     }
 119     diff_buf.create(I0.size(), CV_32FC1);
 120
 121     // create the scales
 122     for (int s = 1; s < nscales; ++s)
 123     {
 124         cuda::resize(I0s[s-1], I0s[s], Size(), scaleStep, scaleStep);
 125         cuda::resize(I1s[s-1], I1s[s], Size(), scaleStep, scaleStep);
 126
 127         if (I0s[s].cols < 16 || I0s[s].rows < 16)
 128         {
 129             nscales = s;
 130             break;
 131         }
 132
 133         if (useInitialFlow)
 134         {
 135             cuda::resize(u1s[s-1], u1s[s], Size(), scaleStep, scaleStep);
 136             cuda::resize(u2s[s-1], u2s[s], Size(), scaleStep, scaleStep);
 137
 138             cuda::multiply(u1s[s], Scalar::all(scaleStep), u1s[s]);
 139             cuda::multiply(u2s[s], Scalar::all(scaleStep), u2s[s]);
 140         }
 141         else
 142         {
 143             u1s[s].create(I0s[s].size(), CV_32FC1);
 144             u2s[s].create(I0s[s].size(), CV_32FC1);
 145         }
 146         if (gamma)
 147             u3s[s].create(I0s[s].size(), CV_32FC1);
 148     }
 149
 150     if (!useInitialFlow)
 151     {
 152         u1s[nscales-1].setTo(Scalar::all(0));
 153         u2s[nscales-1].setTo(Scalar::all(0));
 154     }
 155     if (gamma)
 156         u3s[nscales - 1].setTo(Scalar::all(0));
 157
 158     // pyramidal structure for computing the optical flow
 159     for (int s = nscales - 1; s >= 0; --s)
 160     {
 161         // compute the optical flow at the current scale
 162         procOneScale(I0s[s], I1s[s], u1s[s], u2s[s], u3s[s]);
 163
 164         // if this was the last scale, finish now
 165         if (s == 0)
 166             break;
 167
 168         // otherwise, upsample the optical flow
 169
 170         // zoom the optical flow for the next finer scale
 171         cuda::resize(u1s[s], u1s[s - 1], I0s[s - 1].size());
 172         cuda::resize(u2s[s], u2s[s - 1], I0s[s - 1].size());
 173         if (gamma)
 174             cuda::resize(u3s[s], u3s[s - 1], I0s[s - 1].size());
 175
 176         // scale the optical flow with the appropriate zoom factor
 177         cuda::multiply(u1s[s - 1], Scalar::all(1/scaleStep), u1s[s - 1]);
 178         cuda::multiply(u2s[s - 1], Scalar::all(1/scaleStep), u2s[s - 1]);
 179     }
 180 }
 181
 182 namespace tvl1flow
 183 {
 184     void centeredGradient(PtrStepSzf src, PtrStepSzf dx, PtrStepSzf dy);
 185     void warpBackward(PtrStepSzf I0, PtrStepSzf I1, PtrStepSzf I1x, PtrStepSzf I1y, PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf I1w, PtrStepSzf I1wx, PtrStepSzf I1wy, PtrStepSzf grad, PtrStepSzf rho);
 186     void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
 187                    PtrStepSzf grad, PtrStepSzf rho_c,
 188                    PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
 189                    PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf error,
 190                    float l_t, float theta, float gamma, bool calcError);
 191     void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32, float taut, const float gamma);
 192 }
 193
 194 void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2, GpuMat& u3)
 195 {
 196     using namespace tvl1flow;
 197
 198     const double scaledEpsilon = epsilon * epsilon * I0.size().area();
 199
 200     CV_DbgAssert( I1.size() == I0.size() );
 201     CV_DbgAssert( I1.type() == I0.type() );
 202     CV_DbgAssert( u1.size() == I0.size() );
 203     CV_DbgAssert( u2.size() == u1.size() );
 204
 205     GpuMat I1x = I1x_buf(Rect(0, 0, I0.cols, I0.rows));
 206     GpuMat I1y = I1y_buf(Rect(0, 0, I0.cols, I0.rows));
 207     centeredGradient(I1, I1x, I1y);
 208
 209     GpuMat I1w = I1w_buf(Rect(0, 0, I0.cols, I0.rows));
 210     GpuMat I1wx = I1wx_buf(Rect(0, 0, I0.cols, I0.rows));
 211     GpuMat I1wy = I1wy_buf(Rect(0, 0, I0.cols, I0.rows));
 212
 213     GpuMat grad = grad_buf(Rect(0, 0, I0.cols, I0.rows));
 214     GpuMat rho_c = rho_c_buf(Rect(0, 0, I0.cols, I0.rows));
 215
 216     GpuMat p11 = p11_buf(Rect(0, 0, I0.cols, I0.rows));
 217     GpuMat p12 = p12_buf(Rect(0, 0, I0.cols, I0.rows));
 218     GpuMat p21 = p21_buf(Rect(0, 0, I0.cols, I0.rows));
 219     GpuMat p22 = p22_buf(Rect(0, 0, I0.cols, I0.rows));
 220     GpuMat p31, p32;
 221     if (gamma)
 222     {
 223         p31 = p31_buf(Rect(0, 0, I0.cols, I0.rows));
 224         p32 = p32_buf(Rect(0, 0, I0.cols, I0.rows));
 225     }
 226     p11.setTo(Scalar::all(0));
 227     p12.setTo(Scalar::all(0));
 228     p21.setTo(Scalar::all(0));
 229     p22.setTo(Scalar::all(0));
 230     if (gamma)
 231     {
 232         p31.setTo(Scalar::all(0));
 233         p32.setTo(Scalar::all(0));
 234     }
 235
 236     GpuMat diff = diff_buf(Rect(0, 0, I0.cols, I0.rows));
 237
 238     const float l_t = static_cast<float>(lambda * theta);
 239     const float taut = static_cast<float>(tau / theta);
 240
 241     for (int warpings = 0; warpings < warps; ++warpings)
 242     {
 243         warpBackward(I0, I1, I1x, I1y, u1, u2, I1w, I1wx, I1wy, grad, rho_c);
 244
 245         double error = std::numeric_limits<double>::max();
 246         double prevError = 0.0;
 247         for (int n = 0; error > scaledEpsilon && n < iterations; ++n)
 248         {
 249             // some tweaks to make sum operation less frequently
 250             bool calcError = (epsilon > 0) && (n & 0x1) && (prevError < scaledEpsilon);
 251             cv::Mat m1(u3);
 252             estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, p31, p32, u1, u2, u3, diff, l_t, static_cast<float>(theta), gamma, calcError);
 253             if (calcError)
 254             {
 255                 error = cuda::sum(diff, norm_buf)[0];
 256                 prevError = error;
 257             }
 258             else
 259             {
 260                 error = std::numeric_limits<double>::max();
 261                 prevError -= scaledEpsilon;
 262             }
 263
 264             estimateDualVariables(u1, u2, u3, p11, p12, p21, p22, p31, p32, taut, gamma);
 265         }
 266     }
 267 }
 268
 269 void cv::cuda::OpticalFlowDual_TVL1_CUDA::collectGarbage()
 270 {
 271     I0s.clear();
 272     I1s.clear();
 273     u1s.clear();
 274     u2s.clear();
 275     u3s.clear();
 276
 277     I1x_buf.release();
 278     I1y_buf.release();
 279
 280     I1w_buf.release();
 281     I1wx_buf.release();
 282     I1wy_buf.release();
 283
 284     grad_buf.release();
 285     rho_c_buf.release();
 286
 287     p11_buf.release();
 288     p12_buf.release();
 289     p21_buf.release();
 290     p22_buf.release();
 291     if (gamma)
 292     {
 293         p31_buf.release();
 294         p32_buf.release();
 295     }
 296     diff_buf.release();
 297     norm_buf.release();
 298 }
 299
 300 #endif // !defined HAVE_CUDA || defined(CUDA_DISABLER)