Added linear programming based stabilizer (videostab)
authorAlexey Spizhevoy <no@email>
Mon, 16 Apr 2012 06:41:06 +0000 (06:41 +0000)
committerAlexey Spizhevoy <no@email>
Mon, 16 Apr 2012 06:41:06 +0000 (06:41 +0000)
modules/videostab/include/opencv2/videostab/global_motion.hpp
modules/videostab/include/opencv2/videostab/inpainting.hpp
modules/videostab/include/opencv2/videostab/motion_stabilizing.hpp
modules/videostab/src/global_motion.cpp
modules/videostab/src/inpainting.cpp
modules/videostab/src/motion_stabilizing.cpp
modules/videostab/src/stabilizer.cpp
modules/videostab/src/wobble_suppression.cpp
samples/cpp/videostab.cpp

index c523c5a..3f9b3d0 100644 (file)
@@ -57,16 +57,16 @@ namespace videostab
 
 enum MotionModel
 {
-    TRANSLATION = 0,
-    TRANSLATION_AND_SCALE = 1,
-    LINEAR_SIMILARITY = 2,
-    AFFINE = 3,
-    HOMOGRAPHY = 4,
-    UNKNOWN = 5
+    MM_TRANSLATION = 0,
+    MM_TRANSLATION_AND_SCALE = 1,
+    MM_LINEAR_SIMILARITY = 2,
+    MM_AFFINE = 3,
+    MM_HOMOGRAPHY = 4,
+    MM_UNKNOWN = 5
 };
 
 CV_EXPORTS Mat estimateGlobalMotionLeastSquares(
-        int npoints, Point2f *points0, Point2f *points1, int model = AFFINE, float *rmse = 0);
+        int npoints, Point2f *points0, Point2f *points1, int model = MM_AFFINE, float *rmse = 0);
 
 struct CV_EXPORTS RansacParams
 {
@@ -81,14 +81,14 @@ struct CV_EXPORTS RansacParams
 
     static RansacParams default2dMotion(MotionModel model)
     {
-        CV_Assert(model < UNKNOWN);
-        if (model == TRANSLATION)
+        CV_Assert(model < MM_UNKNOWN);
+        if (model == MM_TRANSLATION)
             return RansacParams(1, 0.5f, 0.5f, 0.99f);
-        if (model == TRANSLATION_AND_SCALE)
+        if (model == MM_TRANSLATION_AND_SCALE)
             return RansacParams(2, 0.5f, 0.5f, 0.99f);
-        if (model == LINEAR_SIMILARITY)
+        if (model == MM_LINEAR_SIMILARITY)
             return RansacParams(2, 0.5f, 0.5f, 0.99f);
-        if (model == AFFINE)
+        if (model == MM_AFFINE)
             return RansacParams(3, 0.5f, 0.5f, 0.99f);
         return RansacParams(4, 0.5f, 0.5f, 0.99f);
     }
@@ -96,13 +96,13 @@ struct CV_EXPORTS RansacParams
 
 CV_EXPORTS Mat estimateGlobalMotionRobust(
         const std::vector<Point2f> &points0, const std::vector<Point2f> &points1,
-        int model = AFFINE, const RansacParams &params = RansacParams::default2dMotion(AFFINE),
+        int model = MM_AFFINE, const RansacParams &params = RansacParams::default2dMotion(MM_AFFINE),
         float *rmse = 0, int *ninliers = 0);
 
 class CV_EXPORTS GlobalMotionEstimatorBase
 {
 public:
-    GlobalMotionEstimatorBase() : motionModel_(UNKNOWN) {}
+    GlobalMotionEstimatorBase() : motionModel_(MM_UNKNOWN) {}
     virtual ~GlobalMotionEstimatorBase() {}
 
     virtual void setMotionModel(MotionModel val) { motionModel_ = val; }
@@ -138,7 +138,7 @@ private:
 class CV_EXPORTS PyrLkRobustMotionEstimator : public GlobalMotionEstimatorBase
 {
 public:
-    PyrLkRobustMotionEstimator(MotionModel model = AFFINE);
+    PyrLkRobustMotionEstimator(MotionModel model = MM_AFFINE);
 
     void setDetector(Ptr<FeatureDetector> val) { detector_ = val; }
     Ptr<FeatureDetector> detector() const { return detector_; }
index 193c295..580413a 100644 (file)
@@ -59,7 +59,7 @@ class CV_EXPORTS InpainterBase
 {
 public:
     InpainterBase()
-        : radius_(0), motionModel_(UNKNOWN), frames_(0), motions_(0),
+        : radius_(0), motionModel_(MM_UNKNOWN), frames_(0), motions_(0),
           stabilizedFrames_(0), stabilizationMotions_(0) {}
 
     virtual ~InpainterBase() {}
index 6d1522b..78c9f52 100644 (file)
@@ -61,7 +61,7 @@ public:
     // assumes that [0, size-1) is in or equals to [range.first, range.second)
     virtual void stabilize(
             int size, const std::vector<Mat> &motions, std::pair<int,int> range,
-            Mat *stabilizationMotions) const = 0;
+            Mat *stabilizationMotions) = 0;
 };
 
 class CV_EXPORTS MotionStabilizationPipeline : public IMotionStabilizer
@@ -72,7 +72,7 @@ public:
 
     virtual void stabilize(
             int size, const std::vector<Mat> &motions, std::pair<int,int> range,
-            Mat *stabilizationMotions) const;
+            Mat *stabilizationMotions);
 
 private:
     std::vector<Ptr<IMotionStabilizer> > stabilizers_;
@@ -84,11 +84,11 @@ public:
     virtual ~MotionFilterBase() {}
 
     virtual Mat stabilize(
-            int idx, const std::vector<Mat> &motions, std::pair<int,int> range) const = 0;
+            int idx, const std::vector<Mat> &motions, std::pair<int,int> range) = 0;
 
     virtual void stabilize(
             int size, const std::vector<Mat> &motions, std::pair<int,int> range,
-            Mat *stabilizationMotions) const;
+            Mat *stabilizationMotions);
 };
 
 class CV_EXPORTS GaussianMotionFilter : public MotionFilterBase
@@ -101,7 +101,7 @@ public:
     float stdev() const { return stdev_; }
 
     virtual Mat stabilize(
-            int idx, const std::vector<Mat> &motions, std::pair<int,int> range) const;
+            int idx, const std::vector<Mat> &motions, std::pair<int,int> range);
 
 private:
     int radius_;
@@ -112,17 +112,49 @@ private:
 class CV_EXPORTS LpMotionStabilizer : public IMotionStabilizer
 {
 public:
-    LpMotionStabilizer(MotionModel model = LINEAR_SIMILARITY);
+    LpMotionStabilizer(MotionModel model = MM_LINEAR_SIMILARITY);
 
     void setMotionModel(MotionModel val) { model_ = val; }
     MotionModel motionModel() const { return model_; }
 
+    void setFrameSize(Size val) { frameSize_ = val; }
+    Size frameSize() const { return frameSize_; }
+
+    void setTrimRatio(float val) { trimRatio_ = val; }
+    float trimRatio() const { return trimRatio_; }
+
+    void setWeight1(float val) { w1_ = val; }
+    float weight1() const { return w1_; }
+
+    void setWeight2(float val) { w2_ = val; }
+    float weight2() const { return w2_; }
+
+    void setWeight3(float val) { w3_ = val; }
+    float weight3() const { return w3_; }
+
+    void setWeight4(float val) { w4_ = val; }
+    float weight4() const { return w4_; }
+
     virtual void stabilize(
             int size, const std::vector<Mat> &motions, std::pair<int,int> range,
-            Mat *stabilizationMotions) const;
+            Mat *stabilizationMotions);
 
 private:
     MotionModel model_;
+    Size frameSize_;
+    float trimRatio_;
+    float w1_, w2_, w3_, w4_;
+
+    std::vector<double> obj_, collb_, colub_;
+    std::vector<int> rows_, cols_;
+    std::vector<double> elems_, rowlb_, rowub_;
+
+    void set(int row, int col, double coef)
+    {
+        rows_.push_back(row);
+        cols_.push_back(col);
+        elems_.push_back(coef);
+    }
 };
 
 CV_EXPORTS Mat ensureInclusionConstraint(const Mat &M, Size size, float trimRatio);
index c419324..806abad 100644 (file)
@@ -228,7 +228,7 @@ static Mat estimateGlobMotionLeastSquaresAffine(
 Mat estimateGlobalMotionLeastSquares(
         int npoints, Point2f *points0, Point2f *points1, int model, float *rmse)
 {
-    CV_Assert(model <= AFFINE);
+    CV_Assert(model <= MM_AFFINE);
 
     typedef Mat (*Impl)(int, Point2f*, Point2f*, float*);
     static Impl impls[] = { estimateGlobMotionLeastSquaresTranslation,
@@ -244,7 +244,7 @@ Mat estimateGlobalMotionRobust(
         const vector<Point2f> &points0, const vector<Point2f> &points1, int model,
         const RansacParams &params, float *rmse, int *ninliers)
 {
-    CV_Assert(model <= AFFINE);
+    CV_Assert(model <= MM_AFFINE);
     CV_Assert(points0.size() == points1.size());
 
     const int npoints = static_cast<int>(points0.size());
@@ -436,7 +436,7 @@ Mat PyrLkRobustMotionEstimator::estimate(const Mat &frame0, const Mat &frame1, b
     int ninliers;
     Mat_<float> M;
 
-    if (motionModel_ != HOMOGRAPHY)
+    if (motionModel_ != MM_HOMOGRAPHY)
         M = estimateGlobalMotionRobust(
                 pointsPrevGood_, pointsGood_, motionModel_, ransacParams_, &rmse, &ninliers);
     else
index ab7348a..9eed078 100644 (file)
@@ -375,7 +375,7 @@ void MotionInpainter::inpaint(int idx, Mat &frame, Mat &mask)
 
         frame1_ = at(neighbor, *frames_);        
 
-        if (motionModel_ != HOMOGRAPHY)
+        if (motionModel_ != MM_HOMOGRAPHY)
             warpAffine(
                     frame1_, transformedFrame1_, motion1to0(Rect(0,0,3,2)), frame1_.size(),
                     INTER_LINEAR, borderMode_);
@@ -388,7 +388,7 @@ void MotionInpainter::inpaint(int idx, Mat &frame, Mat &mask)
 
         // warp mask
 
-        if (motionModel_ != HOMOGRAPHY)
+        if (motionModel_ != MM_HOMOGRAPHY)
             warpAffine(
                     mask1_, transformedMask1_, motion1to0(Rect(0,0,3,2)), mask1_.size(),
                     INTER_NEAREST);
index e0cf55a..27a4595 100644 (file)
 
 #ifdef HAVE_CLP
 #include "coin/ClpSimplex.hpp"
+#include "coin/ClpPresolve.hpp"
+#include "coin/ClpPrimalColumnSteepest.hpp"
+#include "coin/ClpDualRowSteepest.hpp"
 #endif
 
+#define INF 1e10
+
 using namespace std;
 
 namespace cv
@@ -58,7 +63,7 @@ namespace videostab
 
 void MotionStabilizationPipeline::stabilize(
         int size, const vector<Mat> &motions, pair<int,int> range,
-        Mat *stabilizationMotions) const
+        Mat *stabilizationMotions)
 {
     vector<Mat> updatedMotions(motions.size());
     for (size_t i = 0; i < motions.size(); ++i)
@@ -87,7 +92,7 @@ void MotionStabilizationPipeline::stabilize(
 
 
 void MotionFilterBase::stabilize(
-        int size, const vector<Mat> &motions, pair<int,int> range, Mat *stabilizationMotions) const
+        int size, const vector<Mat> &motions, pair<int,int> range, Mat *stabilizationMotions)
 {
     for (int i = 0; i < size; ++i)
         stabilizationMotions[i] = stabilize(i, motions, range);
@@ -108,7 +113,7 @@ void GaussianMotionFilter::setParams(int radius, float stdev)
 }
 
 
-Mat GaussianMotionFilter::stabilize(int idx, const vector<Mat> &motions, pair<int,int> range) const
+Mat GaussianMotionFilter::stabilize(int idx, const vector<Mat> &motions, pair<int,int> range)
 {
     const Mat &cur = at(idx, motions);
     Mat res = Mat::zeros(cur.size(), cur.type());
@@ -260,23 +265,449 @@ float estimateOptimalTrimRatio(const Mat &M, Size size)
 LpMotionStabilizer::LpMotionStabilizer(MotionModel model)
 {
     setMotionModel(model);
+    setFrameSize(Size(0,0));
+    setTrimRatio(0.1);
+    setWeight1(1);
+    setWeight2(10);
+    setWeight3(100);
+    setWeight4(100);
 }
 
 
 #ifndef HAVE_CLP
+
 void LpMotionStabilizer::stabilize(int, const vector<Mat>&, pair<int,int>, Mat*) const
 {
     CV_Error(CV_StsError, "The library is built without Clp support");
 }
+
 #else
+
 void LpMotionStabilizer::stabilize(
         int size, const vector<Mat> &motions, pair<int,int> range,
-        Mat *stabilizationMotions) const
+        Mat *stabilizationMotions)
 {
-    // TODO implement
-    CV_Error(CV_StsNotImplemented, "LpMotionStabilizer::stabilize");
+    CV_Assert(model_ == MM_LINEAR_SIMILARITY);
+
+    int N = size;
+    const vector<Mat> &M = motions;
+    Mat *S = stabilizationMotions;
+
+    double w = frameSize_.width, h = frameSize_.height;
+    double tw = w * trimRatio_, th = h * trimRatio_;
+
+    int ncols = 4*N + 6*(N-1) + 6*(N-2) + 6*(N-3);
+    int nrows = 8*N + 2*6*(N-1) + 2*6*(N-2) + 2*6*(N-3);
+
+    obj_.assign(ncols, 0);
+    collb_.assign(ncols, -INF);
+    colub_.assign(ncols, INF);
+    int c = 4*N;
+
+    // for each slack variable e[t] (error bound)
+    for (int t = 0; t < N-1; ++t, c += 6)
+    {
+        // e[t](0,0)
+        obj_[c] = w4_*w1_;
+        collb_[c] = 0;
+
+        // e[t](0,1)
+        obj_[c+1] = w4_*w1_;
+        collb_[c+1] = 0;
+
+        // e[t](0,2)
+        obj_[c+2] = w1_;
+        collb_[c+2] = 0;
+
+        // e[t](1,0)
+        obj_[c+3] = w4_*w1_;
+        collb_[c+3] = 0;
+
+        // e[t](1,1)
+        obj_[c+4] = w4_*w1_;
+        collb_[c+4] = 0;
+
+        // e[t](1,2)
+        obj_[c+5] = w1_;
+        collb_[c+5] = 0;
+    }
+    for (int t = 0; t < N-2; ++t, c += 6)
+    {
+        // e[t](0,0)
+        obj_[c] = w4_*w2_;
+        collb_[c] = 0;
+
+        // e[t](0,1)
+        obj_[c+1] = w4_*w2_;
+        collb_[c+1] = 0;
+
+        // e[t](0,2)
+        obj_[c+2] = w2_;
+        collb_[c+2] = 0;
+
+        // e[t](1,0)
+        obj_[c+3] = w4_*w2_;
+        collb_[c+3] = 0;
+
+        // e[t](1,1)
+        obj_[c+4] = w4_*w2_;
+        collb_[c+4] = 0;
+
+        // e[t](1,2)
+        obj_[c+5] = w2_;
+        collb_[c+5] = 0;
+    }
+    for (int t = 0; t < N-3; ++t, c += 6)
+    {
+        // e[t](0,0)
+        obj_[c] = w4_*w3_;
+        collb_[c] = 0;
+
+        // e[t](0,1)
+        obj_[c+1] = w4_*w3_;
+        collb_[c+1] = 0;
+
+        // e[t](0,2)
+        obj_[c+2] = w3_;
+        collb_[c+2] = 0;
+
+        // e[t](1,0)
+        obj_[c+3] = w4_*w3_;
+        collb_[c+3] = 0;
+
+        // e[t](1,1)
+        obj_[c+4] = w4_*w3_;
+        collb_[c+4] = 0;
+
+        // e[t](1,2)
+        obj_[c+5] = w3_;
+        collb_[c+5] = 0;
+    }
+
+    elems_.clear();
+    rowlb_.assign(nrows, -INF);
+    rowub_.assign(nrows, INF);
+
+    vector<CoinShallowPackedVector> packedRows;
+    packedRows.reserve(nrows);
+
+    int r = 0;
+
+    // frame corners
+    const Point2d pt[] = {Point2d(0,0), Point2d(w,0), Point2d(w,h), Point2d(0,h)};
+
+    // for each frame
+    for (int t = 0; t < N; ++t)
+    {
+        c = 4*t;
+
+        // for each frame corner
+        for (int i = 0; i < 4; ++i, r += 2)
+        {
+            set(r, c, pt[i].x); set(r, c+1, pt[i].y); set(r, c+2, 1);
+            set(r+1, c, pt[i].y); set(r+1, c+1, -pt[i].x); set(r+1, c+3, 1);
+            rowlb_[r] = pt[i].x-tw; rowub_[r] = pt[i].x+tw;
+            rowlb_[r+1] = pt[i].y-th; rowub_[r+1] = pt[i].y+th;
+        }
+    }
+
+    // for each S[t+1]M[t] - S[t] - e[t] <= 0 condition
+    for (int t = 0; t < N-1; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M);
+
+        c = 4*t;
+        set(r, c, -1);
+        set(r+1, c+1, -1);
+        set(r+2, c+2, -1);
+        set(r+3, c+1, 1);
+        set(r+4, c, -1);
+        set(r+5, c+3, -1);
+
+        c = 4*(t+1);
+        set(r, c, M0(0,0)); set(r, c+1, M0(1,0));
+        set(r+1, c, M0(0,1)); set(r+1, c+1, M0(1,1));
+        set(r+2, c, M0(0,2)); set(r+2, c+1, M0(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M0(1,0)); set(r+3, c+1, -M0(0,0));
+        set(r+4, c, M0(1,1)); set(r+4, c+1, -M0(0,1));
+        set(r+5, c, M0(1,2)); set(r+5, c+1, -M0(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, -1);
+
+        rowub_[r] = 0;
+        rowub_[r+1] = 0;
+        rowub_[r+2] = 0;
+        rowub_[r+3] = 0;
+        rowub_[r+4] = 0;
+        rowub_[r+5] = 0;
+    }
+
+    // for each 0 <= S[t+1]M[t] - S[t] + e[t] condition
+    for (int t = 0; t < N-1; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M);
+
+        c = 4*t;
+        set(r, c, -1);
+        set(r+1, c+1, -1);
+        set(r+2, c+2, -1);
+        set(r+3, c+1, 1);
+        set(r+4, c, -1);
+        set(r+5, c+3, -1);
+
+        c = 4*(t+1);
+        set(r, c, M0(0,0)); set(r, c+1, M0(1,0));
+        set(r+1, c, M0(0,1)); set(r+1, c+1, M0(1,1));
+        set(r+2, c, M0(0,2)); set(r+2, c+1, M0(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M0(1,0)); set(r+3, c+1, -M0(0,0));
+        set(r+4, c, M0(1,1)); set(r+4, c+1, -M0(0,1));
+        set(r+5, c, M0(1,2)); set(r+5, c+1, -M0(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, 1);
+
+        rowlb_[r] = 0;
+        rowlb_[r+1] = 0;
+        rowlb_[r+2] = 0;
+        rowlb_[r+3] = 0;
+        rowlb_[r+4] = 0;
+        rowlb_[r+5] = 0;
+    }
+
+    // for each S[t+2]M[t+1] - S[t+1]*(I+M[t]) + S[t] - e[t] <= 0 condition
+    for (int t = 0; t < N-2; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M), M1 = at(t+1,M);
+
+        c = 4*t;
+        set(r, c, 1);
+        set(r+1, c+1, 1);
+        set(r+2, c+2, 1);
+        set(r+3, c+1, -1);
+        set(r+4, c, 1);
+        set(r+5, c+3, 1);
+
+        c = 4*(t+1);
+        set(r, c, -M0(0,0)-1); set(r, c+1, -M0(1,0));
+        set(r+1, c, -M0(0,1)); set(r+1, c+1, -M0(1,1)-1);
+        set(r+2, c, -M0(0,2)); set(r+2, c+1, -M0(1,2)); set(r+2, c+2, -2);
+        set(r+3, c, -M0(1,0)); set(r+3, c+1, M0(0,0)+1);
+        set(r+4, c, -M0(1,1)-1); set(r+4, c+1, M0(0,1));
+        set(r+5, c, -M0(1,2)); set(r+5, c+1, M0(0,2)); set(r+5, c+3, -2);
+
+        c = 4*(t+2);
+        set(r, c, M1(0,0)); set(r, c+1, M1(1,0));
+        set(r+1, c, M1(0,1)); set(r+1, c+1, M1(1,1));
+        set(r+2, c, M1(0,2)); set(r+2, c+1, M1(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M1(1,0)); set(r+3, c+1, -M1(0,0));
+        set(r+4, c, M1(1,1)); set(r+4, c+1, -M1(0,1));
+        set(r+5, c, M1(1,2)); set(r+5, c+1, -M1(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*(N-1) + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, -1);
+
+        rowub_[r] = 0;
+        rowub_[r+1] = 0;
+        rowub_[r+2] = 0;
+        rowub_[r+3] = 0;
+        rowub_[r+4] = 0;
+        rowub_[r+5] = 0;
+    }
+
+    // for each 0 <= S[t+2]M[t+1]] - S[t+1]*(I+M[t]) + S[t] + e[t] condition
+    for (int t = 0; t < N-2; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M), M1 = at(t+1,M);
+
+        c = 4*t;
+        set(r, c, 1);
+        set(r+1, c+1, 1);
+        set(r+2, c+2, 1);
+        set(r+3, c+1, -1);
+        set(r+4, c, 1);
+        set(r+5, c+3, 1);
+
+        c = 4*(t+1);
+        set(r, c, -M0(0,0)-1); set(r, c+1, -M0(1,0));
+        set(r+1, c, -M0(0,1)); set(r+1, c+1, -M0(1,1)-1);
+        set(r+2, c, -M0(0,2)); set(r+2, c+1, -M0(1,2)); set(r+2, c+2, -2);
+        set(r+3, c, -M0(1,0)); set(r+3, c+1, M0(0,0)+1);
+        set(r+4, c, -M0(1,1)-1); set(r+4, c+1, M0(0,1));
+        set(r+5, c, -M0(1,2)); set(r+5, c+1, M0(0,2)); set(r+5, c+3, -2);
+
+        c = 4*(t+2);
+        set(r, c, M1(0,0)); set(r, c+1, M1(1,0));
+        set(r+1, c, M1(0,1)); set(r+1, c+1, M1(1,1));
+        set(r+2, c, M1(0,2)); set(r+2, c+1, M1(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M1(1,0)); set(r+3, c+1, -M1(0,0));
+        set(r+4, c, M1(1,1)); set(r+4, c+1, -M1(0,1));
+        set(r+5, c, M1(1,2)); set(r+5, c+1, -M1(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*(N-1) + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, 1);
+
+        rowlb_[r] = 0;
+        rowlb_[r+1] = 0;
+        rowlb_[r+2] = 0;
+        rowlb_[r+3] = 0;
+        rowlb_[r+4] = 0;
+        rowlb_[r+5] = 0;
+    }
+
+    // for each S[t+3]M[t+2] - S[t+2]*(I+2M[t+1]) + S[t+1]*(2*I+M[t]) - S[t] - e[t] <= 0 condition
+    for (int t = 0; t < N-3; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M), M1 = at(t+1,M), M2 = at(t+2,M);
+
+        c = 4*t;
+        set(r, c, -1);
+        set(r+1, c+1, -1);
+        set(r+2, c+2, -1);
+        set(r+3, c+1, 1);
+        set(r+4, c, -1);
+        set(r+5, c+3, -1);
+
+        c = 4*(t+1);
+        set(r, c, M0(0,0)+2); set(r, c+1, M0(1,0));
+        set(r+1, c, M0(0,1)); set(r+1, c+1, M0(1,1)+2);
+        set(r+2, c, M0(0,2)); set(r+2, c+1, M0(1,2)); set(r+2, c+2, 3);
+        set(r+3, c, M0(1,0)); set(r+3, c+1, -M0(0,0)-2);
+        set(r+4, c, M0(1,1)+2); set(r+4, c+1, -M0(0,1));
+        set(r+5, c, M0(1,2)); set(r+5, c+1, -M0(0,2)); set(r+5, c+3, 3);
+
+        c = 4*(t+2);
+        set(r, c, -2*M1(0,0)-1); set(r, c+1, -2*M1(1,0));
+        set(r+1, c, -2*M1(0,1)); set(r+1, c+1, -2*M1(1,1)-1);
+        set(r+2, c, -2*M1(0,2)); set(r+2, c+1, -2*M1(1,2)); set(r+2, c+2, -3);
+        set(r+3, c, -2*M1(1,0)); set(r+3, c+1, 2*M1(0,0)+1);
+        set(r+4, c, -2*M1(1,1)-1); set(r+4, c+1, 2*M1(0,1));
+        set(r+5, c, -2*M1(1,2)); set(r+5, c+1, 2*M1(0,2)); set(r+5, c+3, -3);
+
+        c = 4*(t+3);
+        set(r, c, M2(0,0)); set(r, c+1, M2(1,0));
+        set(r+1, c, M2(0,1)); set(r+1, c+1, M2(1,1));
+        set(r+2, c, M2(0,2)); set(r+2, c+1, M2(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M2(1,0)); set(r+3, c+1, -M2(0,0));
+        set(r+4, c, M2(1,1)); set(r+4, c+1, -M2(0,1));
+        set(r+5, c, M2(1,2)); set(r+5, c+1, -M2(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*(N-1) + 6*(N-2) + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, -1);
+
+        rowub_[r] = 0;
+        rowub_[r+1] = 0;
+        rowub_[r+2] = 0;
+        rowub_[r+3] = 0;
+        rowub_[r+4] = 0;
+        rowub_[r+5] = 0;
+    }
+
+    // for each 0 <= S[t+3]M[t+2] - S[t+2]*(I+2M[t+1]) + S[t+1]*(2*I+M[t]) + e[t] condition
+    for (int t = 0; t < N-3; ++t, r += 6)
+    {
+        Mat_<float> M0 = at(t,M), M1 = at(t+1,M), M2 = at(t+2,M);
+
+        c = 4*t;
+        set(r, c, -1);
+        set(r+1, c+1, -1);
+        set(r+2, c+2, -1);
+        set(r+3, c+1, 1);
+        set(r+4, c, -1);
+        set(r+5, c+3, -1);
+
+        c = 4*(t+1);
+        set(r, c, M0(0,0)+2); set(r, c+1, M0(1,0));
+        set(r+1, c, M0(0,1)); set(r+1, c+1, M0(1,1)+2);
+        set(r+2, c, M0(0,2)); set(r+2, c+1, M0(1,2)); set(r+2, c+2, 3);
+        set(r+3, c, M0(1,0)); set(r+3, c+1, -M0(0,0)-2);
+        set(r+4, c, M0(1,1)+2); set(r+4, c+1, -M0(0,1));
+        set(r+5, c, M0(1,2)); set(r+5, c+1, -M0(0,2)); set(r+5, c+3, 3);
+
+        c = 4*(t+2);
+        set(r, c, -2*M1(0,0)-1); set(r, c+1, -2*M1(1,0));
+        set(r+1, c, -2*M1(0,1)); set(r+1, c+1, -2*M1(1,1)-1);
+        set(r+2, c, -2*M1(0,2)); set(r+2, c+1, -2*M1(1,2)); set(r+2, c+2, -3);
+        set(r+3, c, -2*M1(1,0)); set(r+3, c+1, 2*M1(0,0)+1);
+        set(r+4, c, -2*M1(1,1)-1); set(r+4, c+1, 2*M1(0,1));
+        set(r+5, c, -2*M1(1,2)); set(r+5, c+1, 2*M1(0,2)); set(r+5, c+3, -3);
+
+        c = 4*(t+3);
+        set(r, c, M2(0,0)); set(r, c+1, M2(1,0));
+        set(r+1, c, M2(0,1)); set(r+1, c+1, M2(1,1));
+        set(r+2, c, M2(0,2)); set(r+2, c+1, M2(1,2)); set(r+2, c+2, 1);
+        set(r+3, c, M2(1,0)); set(r+3, c+1, -M2(0,0));
+        set(r+4, c, M2(1,1)); set(r+4, c+1, -M2(0,1));
+        set(r+5, c, M2(1,2)); set(r+5, c+1, -M2(0,2)); set(r+5, c+3, 1);
+
+        c = 4*N + 6*(N-1) + 6*(N-2) + 6*t;
+        for (int i = 0; i < 6; ++i)
+            set(r+i, c+i, 1);
+
+        rowlb_[r] = 0;
+        rowlb_[r+1] = 0;
+        rowlb_[r+2] = 0;
+        rowlb_[r+3] = 0;
+        rowlb_[r+4] = 0;
+        rowlb_[r+5] = 0;
+    }
+
+    // solve
+
+    CoinPackedMatrix A(true, &rows_[0], &cols_[0], &elems_[0], elems_.size());
+    A.setDimensions(nrows, ncols);
+
+    ClpSimplex model;
+    model.loadProblem(A, &collb_[0], &colub_[0], &obj_[0], &rowlb_[0], &rowub_[0]);
+
+    ClpDualRowSteepest dualSteep(1);
+    model.setDualRowPivotAlgorithm(dualSteep);
+
+    ClpPrimalColumnSteepest primalSteep(1);
+    model.setPrimalColumnPivotAlgorithm(primalSteep);
+
+    model.scaling(1);
+
+    ClpPresolve presolveInfo;
+    Ptr<ClpSimplex> presolvedModel = presolveInfo.presolvedModel(model);
+
+    if (!presolvedModel.empty())
+    {
+        presolvedModel->dual();
+        presolveInfo.postsolve(true);
+        model.checkSolution();
+        model.primal(1);
+    }
+    else
+    {
+        model.dual();
+        model.checkSolution();
+        model.primal(1);
+    }
+
+    // save results
+
+    const double *sol = model.getColSolution();
+    c = 0;
+
+    for (int t = 0; t < N; ++t, c += 4)
+    {
+        Mat_<float> S0 = Mat::eye(3, 3, CV_32F);
+        S0(1,1) = S0(0,0) = sol[c];
+        S0(0,1) = sol[c+1];
+        S0(1,0) = -sol[c+1];
+        S0(0,2) = sol[c+2];
+        S0(1,2) = sol[c+3];
+        S[t] = S0;
+    }
 }
-#endif
+
+#endif // #ifndef HAVE_CLP
 
 } // namespace videostab
 } // namespace cv
index 67b7491..57434ab 100644 (file)
@@ -194,7 +194,7 @@ void StabilizerBase::stabilizeFrame()
 
     // apply stabilization transformation
 
-    if (motionEstimator_->motionModel() != HOMOGRAPHY)
+    if (motionEstimator_->motionModel() != MM_HOMOGRAPHY)
         warpAffine(
                 preProcessedFrame_, at(curStabilizedPos_, stabilizedFrames_),
                 stabilizationMotion(Rect(0,0,3,2)), frameSize_, INTER_LINEAR, borderMode_);
@@ -205,7 +205,7 @@ void StabilizerBase::stabilizeFrame()
 
     if (doInpainting_)
     {
-        if (motionEstimator_->motionModel() != HOMOGRAPHY)
+        if (motionEstimator_->motionModel() != MM_HOMOGRAPHY)
             warpAffine(
                     frameMask_, at(curStabilizedPos_, stabilizedMasks_),
                     stabilizationMotion(Rect(0,0,3,2)), frameSize_, INTER_NEAREST);
index 450ec0c..a09437f 100644 (file)
@@ -54,8 +54,8 @@ namespace videostab
 WobbleSuppressorBase::WobbleSuppressorBase() : motions_(0), stabilizationMotions_(0)
 {
     PyrLkRobustMotionEstimator *est = new PyrLkRobustMotionEstimator();
-    est->setMotionModel(HOMOGRAPHY);
-    est->setRansacParams(RansacParams::default2dMotion(HOMOGRAPHY));
+    est->setMotionModel(MM_HOMOGRAPHY);
+    est->setRansacParams(RansacParams::default2dMotion(MM_HOMOGRAPHY));
     setMotionEstimator(est);
 }
 
index f401dcc..dccd413 100644 (file)
@@ -86,6 +86,18 @@ void printHelp()
             "  --stdev=(<float_number>|auto)\n"
             "      Set smoothing weights standard deviation. The default is auto\n"
             "      (i.e. sqrt(radius)).\n"
+            "  -lp, --lp-stab=(yes|no)\n"
+            "      Turn on/off linear programming based stabilization method.\n"
+            "  --lp-trim-ratio=(<float_number>|auto)\n"
+            "      Trimming ratio used in linear programming based method.\n"
+            "  --lp-w1=(<float_number>|1)\n"
+            "      1st derivative weight. The default is 1.\n"
+            "  --lp-w2=(<float_number>|10)\n"
+            "      2nd derivative weight. The default is 10.\n"
+            "  --lp-w3=(<float_number>|100)\n"
+            "      3rd derivative weight. The default is 100.\n"
+            "  --lp-w4=(<float_number>|100)\n"
+            "      Non-translation motion components weight. The default is 100.\n\n"
             "  --deblur=(yes|no)\n"
             "      Do deblurring.\n"
             "  --deblur-sens=<float_number>\n"
@@ -160,6 +172,12 @@ int main(int argc, const char **argv)
                 "{ lm | load-motions | no | }"
                 "{ r | radius | 15 | }"
                 "{ | stdev | auto | }"
+                "{ lp | lp-stab | no | }"
+                "{ | lp-trim-ratio | auto | }"
+                "{ | lp-w1 | 1 | }"
+                "{ | lp-w2 | 10 | }"
+                "{ | lp-w3 | 100 | }"
+                "{ | lp-w4 | 100 | }"
                 "{ | deblur | no | }"
                 "{ | deblur-sens | 0.1 | }"
                 "{ et | est-trim | yes | }"
@@ -194,19 +212,37 @@ int main(int argc, const char **argv)
         {
             printHelp();
             return 0;
-        }               
+        }
+
+        string inputPath = arg("1");
+        if (inputPath.empty()) throw runtime_error("specify video file path");
+
+        VideoFileSource *source = new VideoFileSource(inputPath);
+        cout << "frame count (rough): " << source->count() << endl;
+        if (arg("fps") == "auto") outputFps = source->fps();  else outputFps = argd("fps");
 
         StabilizerBase *stabilizer;
 
         bool isTwoPass =
-                arg("est-trim") == "yes" || arg("wobble-suppress") == "yes";
+                arg("est-trim") == "yes" || arg("wobble-suppress") == "yes" || arg("lp-stab") == "yes";
 
         if (isTwoPass)
         {
             TwoPassStabilizer *twoPassStabilizer = new TwoPassStabilizer();
             stabilizer = twoPassStabilizer;
             twoPassStabilizer->setEstimateTrimRatio(arg("est-trim") == "yes");
-            if (arg("stdev") == "auto")
+            if (arg("lp-stab") == "yes")
+            {
+                LpMotionStabilizer *stab = new LpMotionStabilizer();
+                stab->setFrameSize(Size(source->width(), source->height()));
+                stab->setTrimRatio(arg("lp-trim-ratio") == "auto" ? argf("trim-ratio") : argf("lp-trim-ratio"));
+                stab->setWeight1(argf("lp-w1"));
+                stab->setWeight2(argf("lp-w2"));
+                stab->setWeight3(argf("lp-w3"));
+                stab->setWeight4(argf("lp-w4"));
+                twoPassStabilizer->setMotionStabilizer(stab);
+            }
+            else if (arg("stdev") == "auto")
                 twoPassStabilizer->setMotionStabilizer(new GaussianMotionFilter(argi("radius")));
             else
                 twoPassStabilizer->setMotionStabilizer(new GaussianMotionFilter(argi("radius"), argf("stdev")));
@@ -219,15 +255,15 @@ int main(int argc, const char **argv)
                 PyrLkRobustMotionEstimator *est = 0;
 
                 if (arg("ws-model") == "transl")
-                    est = new PyrLkRobustMotionEstimator(TRANSLATION);
+                    est = new PyrLkRobustMotionEstimator(MM_TRANSLATION);
                 else if (arg("ws-model") == "transl_and_scale")
-                    est = new PyrLkRobustMotionEstimator(TRANSLATION_AND_SCALE);
+                    est = new PyrLkRobustMotionEstimator(MM_TRANSLATION_AND_SCALE);
                 else if (arg("ws-model") == "linear_sim")
-                    est = new PyrLkRobustMotionEstimator(LINEAR_SIMILARITY);
+                    est = new PyrLkRobustMotionEstimator(MM_LINEAR_SIMILARITY);
                 else if (arg("ws-model") == "affine")
-                    est = new PyrLkRobustMotionEstimator(AFFINE);
+                    est = new PyrLkRobustMotionEstimator(MM_AFFINE);
                 else if (arg("ws-model") == "homography")
-                    est = new PyrLkRobustMotionEstimator(HOMOGRAPHY);
+                    est = new PyrLkRobustMotionEstimator(MM_HOMOGRAPHY);
                 else
                 {
                     delete est;
@@ -266,28 +302,22 @@ int main(int argc, const char **argv)
             else
                 onePassStabilizer->setMotionFilter(new GaussianMotionFilter(argi("radius"), argf("stdev")));
         }
-        stabilizedFrames = dynamic_cast<IFrameSource*>(stabilizer);
-
-        string inputPath = arg("1");
-        if (inputPath.empty()) throw runtime_error("specify video file path");
 
-        VideoFileSource *source = new VideoFileSource(inputPath);
-        cout << "frame count (rough): " << source->count() << endl;
-        if (arg("fps") == "auto") outputFps = source->fps();  else outputFps = argd("fps");
         stabilizer->setFrameSource(source);
+        stabilizedFrames = dynamic_cast<IFrameSource*>(stabilizer);
 
         PyrLkRobustMotionEstimator *est = 0;
 
         if (arg("model") == "transl")
-            est = new PyrLkRobustMotionEstimator(TRANSLATION);
+            est = new PyrLkRobustMotionEstimator(MM_TRANSLATION);
         else if (arg("model") == "transl_and_scale")
-            est = new PyrLkRobustMotionEstimator(TRANSLATION_AND_SCALE);
+            est = new PyrLkRobustMotionEstimator(MM_TRANSLATION_AND_SCALE);
         else if (arg("model") == "linear_sim")
-            est = new PyrLkRobustMotionEstimator(LINEAR_SIMILARITY);
+            est = new PyrLkRobustMotionEstimator(MM_LINEAR_SIMILARITY);
         else if (arg("model") == "affine")
-            est = new PyrLkRobustMotionEstimator(AFFINE);
+            est = new PyrLkRobustMotionEstimator(MM_AFFINE);
         else if (arg("model") == "homography")
-            est = new PyrLkRobustMotionEstimator(HOMOGRAPHY);
+            est = new PyrLkRobustMotionEstimator(MM_HOMOGRAPHY);
         else
         {
             delete est;