Home | History | Annotate | Download | only in src
      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 Ptr<OpticalFlowDual_TVL1> cv::cuda::OpticalFlowDual_TVL1::create(double, double, double, int, int, double, int, double, double, bool) { throw_no_cuda(); return Ptr<OpticalFlowDual_TVL1>(); }
     48 
     49 #else
     50 
     51 using namespace cv;
     52 using namespace cv::cuda;
     53 
     54 namespace tvl1flow
     55 {
     56     void centeredGradient(PtrStepSzf src, PtrStepSzf dx, PtrStepSzf dy, cudaStream_t stream);
     57     void warpBackward(PtrStepSzf I0, PtrStepSzf I1, PtrStepSzf I1x, PtrStepSzf I1y,
     58                       PtrStepSzf u1, PtrStepSzf u2,
     59                       PtrStepSzf I1w, PtrStepSzf I1wx, PtrStepSzf I1wy,
     60                       PtrStepSzf grad, PtrStepSzf rho,
     61                       cudaStream_t stream);
     62     void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
     63                    PtrStepSzf grad, PtrStepSzf rho_c,
     64                    PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
     65                    PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3, PtrStepSzf error,
     66                    float l_t, float theta, float gamma, bool calcError,
     67                    cudaStream_t stream);
     68     void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf u3,
     69                                PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, PtrStepSzf p31, PtrStepSzf p32,
     70                                float taut, float gamma,
     71                                cudaStream_t stream);
     72 }
     73 
     74 namespace
     75 {
     76     class OpticalFlowDual_TVL1_Impl : public OpticalFlowDual_TVL1
     77     {
     78     public:
     79         OpticalFlowDual_TVL1_Impl(double tau, double lambda, double theta, int nscales, int warps, double epsilon,
     80                                   int iterations, double scaleStep, double gamma, bool useInitialFlow) :
     81             tau_(tau), lambda_(lambda), gamma_(gamma), theta_(theta), nscales_(nscales), warps_(warps),
     82             epsilon_(epsilon), iterations_(iterations), scaleStep_(scaleStep), useInitialFlow_(useInitialFlow)
     83         {
     84         }
     85 
     86         virtual double getTau() const { return tau_; }
     87         virtual void setTau(double tau) { tau_ = tau; }
     88 
     89         virtual double getLambda() const { return lambda_; }
     90         virtual void setLambda(double lambda) { lambda_ = lambda; }
     91 
     92         virtual double getGamma() const { return gamma_; }
     93         virtual void setGamma(double gamma) { gamma_ = gamma; }
     94 
     95         virtual double getTheta() const { return theta_; }
     96         virtual void setTheta(double theta) { theta_ = theta; }
     97 
     98         virtual int getNumScales() const { return nscales_; }
     99         virtual void setNumScales(int nscales) { nscales_ = nscales; }
    100 
    101         virtual int getNumWarps() const { return warps_; }
    102         virtual void setNumWarps(int warps) { warps_ = warps; }
    103 
    104         virtual double getEpsilon() const { return epsilon_; }
    105         virtual void setEpsilon(double epsilon) { epsilon_ = epsilon; }
    106 
    107         virtual int getNumIterations() const { return iterations_; }
    108         virtual void setNumIterations(int iterations) { iterations_ = iterations; }
    109 
    110         virtual double getScaleStep() const { return scaleStep_; }
    111         virtual void setScaleStep(double scaleStep) { scaleStep_ = scaleStep; }
    112 
    113         virtual bool getUseInitialFlow() const { return useInitialFlow_; }
    114         virtual void setUseInitialFlow(bool useInitialFlow) { useInitialFlow_ = useInitialFlow; }
    115 
    116         virtual void calc(InputArray I0, InputArray I1, InputOutputArray flow, Stream& stream);
    117 
    118     private:
    119         double tau_;
    120         double lambda_;
    121         double gamma_;
    122         double theta_;
    123         int nscales_;
    124         int warps_;
    125         double epsilon_;
    126         int iterations_;
    127         double scaleStep_;
    128         bool useInitialFlow_;
    129 
    130     private:
    131         void calcImpl(const GpuMat& I0, const GpuMat& I1, GpuMat& flowx, GpuMat& flowy, Stream& stream);
    132         void procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2, GpuMat& u3, Stream& stream);
    133 
    134         std::vector<GpuMat> I0s;
    135         std::vector<GpuMat> I1s;
    136         std::vector<GpuMat> u1s;
    137         std::vector<GpuMat> u2s;
    138         std::vector<GpuMat> u3s;
    139 
    140         GpuMat I1x_buf;
    141         GpuMat I1y_buf;
    142 
    143         GpuMat I1w_buf;
    144         GpuMat I1wx_buf;
    145         GpuMat I1wy_buf;
    146 
    147         GpuMat grad_buf;
    148         GpuMat rho_c_buf;
    149 
    150         GpuMat p11_buf;
    151         GpuMat p12_buf;
    152         GpuMat p21_buf;
    153         GpuMat p22_buf;
    154         GpuMat p31_buf;
    155         GpuMat p32_buf;
    156 
    157         GpuMat diff_buf;
    158         GpuMat norm_buf;
    159     };
    160 
    161     void OpticalFlowDual_TVL1_Impl::calc(InputArray _frame0, InputArray _frame1, InputOutputArray _flow, Stream& stream)
    162     {
    163         const GpuMat frame0 = _frame0.getGpuMat();
    164         const GpuMat frame1 = _frame1.getGpuMat();
    165 
    166         BufferPool pool(stream);
    167         GpuMat flowx = pool.getBuffer(frame0.size(), CV_32FC1);
    168         GpuMat flowy = pool.getBuffer(frame0.size(), CV_32FC1);
    169 
    170         calcImpl(frame0, frame1, flowx, flowy, stream);
    171 
    172         GpuMat flows[] = {flowx, flowy};
    173         cuda::merge(flows, 2, _flow, stream);
    174     }
    175 
    176     void OpticalFlowDual_TVL1_Impl::calcImpl(const GpuMat& I0, const GpuMat& I1, GpuMat& flowx, GpuMat& flowy, Stream& stream)
    177     {
    178         CV_Assert( I0.type() == CV_8UC1 || I0.type() == CV_32FC1 );
    179         CV_Assert( I0.size() == I1.size() );
    180         CV_Assert( I0.type() == I1.type() );
    181         CV_Assert( !useInitialFlow_ || (flowx.size() == I0.size() && flowx.type() == CV_32FC1 && flowy.size() == flowx.size() && flowy.type() == flowx.type()) );
    182         CV_Assert( nscales_ > 0 );
    183 
    184         // allocate memory for the pyramid structure
    185         I0s.resize(nscales_);
    186         I1s.resize(nscales_);
    187         u1s.resize(nscales_);
    188         u2s.resize(nscales_);
    189         u3s.resize(nscales_);
    190 
    191         I0.convertTo(I0s[0], CV_32F, I0.depth() == CV_8U ? 1.0 : 255.0, stream);
    192         I1.convertTo(I1s[0], CV_32F, I1.depth() == CV_8U ? 1.0 : 255.0, stream);
    193 
    194         if (!useInitialFlow_)
    195         {
    196             flowx.create(I0.size(), CV_32FC1);
    197             flowy.create(I0.size(), CV_32FC1);
    198         }
    199 
    200         u1s[0] = flowx;
    201         u2s[0] = flowy;
    202         if (gamma_)
    203         {
    204             u3s[0].create(I0.size(), CV_32FC1);
    205         }
    206 
    207         I1x_buf.create(I0.size(), CV_32FC1);
    208         I1y_buf.create(I0.size(), CV_32FC1);
    209 
    210         I1w_buf.create(I0.size(), CV_32FC1);
    211         I1wx_buf.create(I0.size(), CV_32FC1);
    212         I1wy_buf.create(I0.size(), CV_32FC1);
    213 
    214         grad_buf.create(I0.size(), CV_32FC1);
    215         rho_c_buf.create(I0.size(), CV_32FC1);
    216 
    217         p11_buf.create(I0.size(), CV_32FC1);
    218         p12_buf.create(I0.size(), CV_32FC1);
    219         p21_buf.create(I0.size(), CV_32FC1);
    220         p22_buf.create(I0.size(), CV_32FC1);
    221         if (gamma_)
    222         {
    223             p31_buf.create(I0.size(), CV_32FC1);
    224             p32_buf.create(I0.size(), CV_32FC1);
    225         }
    226         diff_buf.create(I0.size(), CV_32FC1);
    227 
    228         // create the scales
    229         for (int s = 1; s < nscales_; ++s)
    230         {
    231             cuda::resize(I0s[s-1], I0s[s], Size(), scaleStep_, scaleStep_, INTER_LINEAR, stream);
    232             cuda::resize(I1s[s-1], I1s[s], Size(), scaleStep_, scaleStep_, INTER_LINEAR, stream);
    233 
    234             if (I0s[s].cols < 16 || I0s[s].rows < 16)
    235             {
    236                 nscales_ = s;
    237                 break;
    238             }
    239 
    240             if (useInitialFlow_)
    241             {
    242                 cuda::resize(u1s[s-1], u1s[s], Size(), scaleStep_, scaleStep_, INTER_LINEAR, stream);
    243                 cuda::resize(u2s[s-1], u2s[s], Size(), scaleStep_, scaleStep_, INTER_LINEAR, stream);
    244 
    245                 cuda::multiply(u1s[s], Scalar::all(scaleStep_), u1s[s], 1, -1, stream);
    246                 cuda::multiply(u2s[s], Scalar::all(scaleStep_), u2s[s], 1, -1, stream);
    247             }
    248             else
    249             {
    250                 u1s[s].create(I0s[s].size(), CV_32FC1);
    251                 u2s[s].create(I0s[s].size(), CV_32FC1);
    252             }
    253             if (gamma_)
    254             {
    255                 u3s[s].create(I0s[s].size(), CV_32FC1);
    256             }
    257         }
    258 
    259         if (!useInitialFlow_)
    260         {
    261             u1s[nscales_-1].setTo(Scalar::all(0), stream);
    262             u2s[nscales_-1].setTo(Scalar::all(0), stream);
    263         }
    264         if (gamma_)
    265         {
    266             u3s[nscales_ - 1].setTo(Scalar::all(0), stream);
    267         }
    268 
    269         // pyramidal structure for computing the optical flow
    270         for (int s = nscales_ - 1; s >= 0; --s)
    271         {
    272             // compute the optical flow at the current scale
    273             procOneScale(I0s[s], I1s[s], u1s[s], u2s[s], u3s[s], stream);
    274 
    275             // if this was the last scale, finish now
    276             if (s == 0)
    277                 break;
    278 
    279             // otherwise, upsample the optical flow
    280 
    281             // zoom the optical flow for the next finer scale
    282             cuda::resize(u1s[s], u1s[s - 1], I0s[s - 1].size(), 0, 0, INTER_LINEAR, stream);
    283             cuda::resize(u2s[s], u2s[s - 1], I0s[s - 1].size(), 0, 0, INTER_LINEAR, stream);
    284             if (gamma_)
    285             {
    286                 cuda::resize(u3s[s], u3s[s - 1], I0s[s - 1].size(), 0, 0, INTER_LINEAR, stream);
    287             }
    288 
    289             // scale the optical flow with the appropriate zoom factor
    290             cuda::multiply(u1s[s - 1], Scalar::all(1/scaleStep_), u1s[s - 1], 1, -1, stream);
    291             cuda::multiply(u2s[s - 1], Scalar::all(1/scaleStep_), u2s[s - 1], 1, -1, stream);
    292         }
    293     }
    294 
    295     void OpticalFlowDual_TVL1_Impl::procOneScale(const GpuMat& I0, const GpuMat& I1, GpuMat& u1, GpuMat& u2, GpuMat& u3, Stream& _stream)
    296     {
    297         using namespace tvl1flow;
    298 
    299         cudaStream_t stream = StreamAccessor::getStream(_stream);
    300 
    301         const double scaledEpsilon = epsilon_ * epsilon_ * I0.size().area();
    302 
    303         CV_DbgAssert( I1.size() == I0.size() );
    304         CV_DbgAssert( I1.type() == I0.type() );
    305         CV_DbgAssert( u1.size() == I0.size() );
    306         CV_DbgAssert( u2.size() == u1.size() );
    307 
    308         GpuMat I1x = I1x_buf(Rect(0, 0, I0.cols, I0.rows));
    309         GpuMat I1y = I1y_buf(Rect(0, 0, I0.cols, I0.rows));
    310         centeredGradient(I1, I1x, I1y, stream);
    311 
    312         GpuMat I1w = I1w_buf(Rect(0, 0, I0.cols, I0.rows));
    313         GpuMat I1wx = I1wx_buf(Rect(0, 0, I0.cols, I0.rows));
    314         GpuMat I1wy = I1wy_buf(Rect(0, 0, I0.cols, I0.rows));
    315 
    316         GpuMat grad = grad_buf(Rect(0, 0, I0.cols, I0.rows));
    317         GpuMat rho_c = rho_c_buf(Rect(0, 0, I0.cols, I0.rows));
    318 
    319         GpuMat p11 = p11_buf(Rect(0, 0, I0.cols, I0.rows));
    320         GpuMat p12 = p12_buf(Rect(0, 0, I0.cols, I0.rows));
    321         GpuMat p21 = p21_buf(Rect(0, 0, I0.cols, I0.rows));
    322         GpuMat p22 = p22_buf(Rect(0, 0, I0.cols, I0.rows));
    323         GpuMat p31, p32;
    324         if (gamma_)
    325         {
    326             p31 = p31_buf(Rect(0, 0, I0.cols, I0.rows));
    327             p32 = p32_buf(Rect(0, 0, I0.cols, I0.rows));
    328         }
    329         p11.setTo(Scalar::all(0), _stream);
    330         p12.setTo(Scalar::all(0), _stream);
    331         p21.setTo(Scalar::all(0), _stream);
    332         p22.setTo(Scalar::all(0), _stream);
    333         if (gamma_)
    334         {
    335             p31.setTo(Scalar::all(0), _stream);
    336             p32.setTo(Scalar::all(0), _stream);
    337         }
    338 
    339         GpuMat diff = diff_buf(Rect(0, 0, I0.cols, I0.rows));
    340 
    341         const float l_t = static_cast<float>(lambda_ * theta_);
    342         const float taut = static_cast<float>(tau_ / theta_);
    343 
    344         for (int warpings = 0; warpings < warps_; ++warpings)
    345         {
    346             warpBackward(I0, I1, I1x, I1y, u1, u2, I1w, I1wx, I1wy, grad, rho_c, stream);
    347 
    348             double error = std::numeric_limits<double>::max();
    349             double prevError = 0.0;
    350             for (int n = 0; error > scaledEpsilon && n < iterations_; ++n)
    351             {
    352                 // some tweaks to make sum operation less frequently
    353                 bool calcError = (epsilon_ > 0) && (n & 0x1) && (prevError < scaledEpsilon);
    354                 estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, p31, p32, u1, u2, u3, diff, l_t, static_cast<float>(theta_), gamma_, calcError, stream);
    355                 if (calcError)
    356                 {
    357                     _stream.waitForCompletion();
    358                     error = cuda::sum(diff, norm_buf)[0];
    359                     prevError = error;
    360                 }
    361                 else
    362                 {
    363                     error = std::numeric_limits<double>::max();
    364                     prevError -= scaledEpsilon;
    365                 }
    366 
    367                 estimateDualVariables(u1, u2, u3, p11, p12, p21, p22, p31, p32, taut, gamma_, stream);
    368             }
    369         }
    370     }
    371 }
    372 
    373 Ptr<OpticalFlowDual_TVL1> cv::cuda::OpticalFlowDual_TVL1::create(
    374             double tau, double lambda, double theta, int nscales, int warps,
    375             double epsilon, int iterations, double scaleStep, double gamma, bool useInitialFlow)
    376 {
    377     return makePtr<OpticalFlowDual_TVL1_Impl>(tau, lambda, theta, nscales, warps,
    378                                               epsilon, iterations, scaleStep, gamma, useInitialFlow);
    379 }
    380 
    381 #endif // !defined HAVE_CUDA || defined(CUDA_DISABLER)
    382