From 40e7990b6bb4847155f56a353101f9270271e314 Mon Sep 17 00:00:00 2001 From: Alexey Spizhevoy Date: Mon, 16 Apr 2012 06:41:06 +0000 Subject: [PATCH] Added linear programming based stabilizer (videostab) --- .../opencv2/videostab/global_motion.hpp | 30 +- .../include/opencv2/videostab/inpainting.hpp | 2 +- .../opencv2/videostab/motion_stabilizing.hpp | 46 +- modules/videostab/src/global_motion.cpp | 6 +- modules/videostab/src/inpainting.cpp | 4 +- modules/videostab/src/motion_stabilizing.cpp | 445 +++++++++++++++++- modules/videostab/src/stabilizer.cpp | 4 +- modules/videostab/src/wobble_suppression.cpp | 4 +- samples/cpp/videostab.cpp | 70 ++- 9 files changed, 552 insertions(+), 59 deletions(-) diff --git a/modules/videostab/include/opencv2/videostab/global_motion.hpp b/modules/videostab/include/opencv2/videostab/global_motion.hpp index c523c5afff..3f9b3d0e85 100644 --- a/modules/videostab/include/opencv2/videostab/global_motion.hpp +++ b/modules/videostab/include/opencv2/videostab/global_motion.hpp @@ -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 &points0, const std::vector &points1, - int model = AFFINE, const RansacParams ¶ms = RansacParams::default2dMotion(AFFINE), + int model = MM_AFFINE, const RansacParams ¶ms = 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 val) { detector_ = val; } Ptr detector() const { return detector_; } diff --git a/modules/videostab/include/opencv2/videostab/inpainting.hpp b/modules/videostab/include/opencv2/videostab/inpainting.hpp index 193c295ddf..580413ad10 100644 --- a/modules/videostab/include/opencv2/videostab/inpainting.hpp +++ b/modules/videostab/include/opencv2/videostab/inpainting.hpp @@ -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() {} diff --git a/modules/videostab/include/opencv2/videostab/motion_stabilizing.hpp b/modules/videostab/include/opencv2/videostab/motion_stabilizing.hpp index 6d1522bbce..78c9f52f58 100644 --- a/modules/videostab/include/opencv2/videostab/motion_stabilizing.hpp +++ b/modules/videostab/include/opencv2/videostab/motion_stabilizing.hpp @@ -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 &motions, std::pair 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 &motions, std::pair range, - Mat *stabilizationMotions) const; + Mat *stabilizationMotions); private: std::vector > stabilizers_; @@ -84,11 +84,11 @@ public: virtual ~MotionFilterBase() {} virtual Mat stabilize( - int idx, const std::vector &motions, std::pair range) const = 0; + int idx, const std::vector &motions, std::pair range) = 0; virtual void stabilize( int size, const std::vector &motions, std::pair 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 &motions, std::pair range) const; + int idx, const std::vector &motions, std::pair 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 &motions, std::pair range, - Mat *stabilizationMotions) const; + Mat *stabilizationMotions); private: MotionModel model_; + Size frameSize_; + float trimRatio_; + float w1_, w2_, w3_, w4_; + + std::vector obj_, collb_, colub_; + std::vector rows_, cols_; + std::vector 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); diff --git a/modules/videostab/src/global_motion.cpp b/modules/videostab/src/global_motion.cpp index c419324ede..806abad7e3 100644 --- a/modules/videostab/src/global_motion.cpp +++ b/modules/videostab/src/global_motion.cpp @@ -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 &points0, const vector &points1, int model, const RansacParams ¶ms, float *rmse, int *ninliers) { - CV_Assert(model <= AFFINE); + CV_Assert(model <= MM_AFFINE); CV_Assert(points0.size() == points1.size()); const int npoints = static_cast(points0.size()); @@ -436,7 +436,7 @@ Mat PyrLkRobustMotionEstimator::estimate(const Mat &frame0, const Mat &frame1, b int ninliers; Mat_ M; - if (motionModel_ != HOMOGRAPHY) + if (motionModel_ != MM_HOMOGRAPHY) M = estimateGlobalMotionRobust( pointsPrevGood_, pointsGood_, motionModel_, ransacParams_, &rmse, &ninliers); else diff --git a/modules/videostab/src/inpainting.cpp b/modules/videostab/src/inpainting.cpp index ab7348a12c..9eed078ba2 100644 --- a/modules/videostab/src/inpainting.cpp +++ b/modules/videostab/src/inpainting.cpp @@ -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); diff --git a/modules/videostab/src/motion_stabilizing.cpp b/modules/videostab/src/motion_stabilizing.cpp index e0cf55aaee..27a45954e3 100644 --- a/modules/videostab/src/motion_stabilizing.cpp +++ b/modules/videostab/src/motion_stabilizing.cpp @@ -47,8 +47,13 @@ #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 &motions, pair range, - Mat *stabilizationMotions) const + Mat *stabilizationMotions) { vector 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 &motions, pair range, Mat *stabilizationMotions) const + int size, const vector &motions, pair 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 &motions, pair range) const +Mat GaussianMotionFilter::stabilize(int idx, const vector &motions, pair 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&, pair, Mat*) const { CV_Error(CV_StsError, "The library is built without Clp support"); } + #else + void LpMotionStabilizer::stabilize( int size, const vector &motions, pair 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 &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 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_ 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_ 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_ 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_ 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_ 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_ 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 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_ 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 diff --git a/modules/videostab/src/stabilizer.cpp b/modules/videostab/src/stabilizer.cpp index 67b749140a..57434ab75f 100644 --- a/modules/videostab/src/stabilizer.cpp +++ b/modules/videostab/src/stabilizer.cpp @@ -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); diff --git a/modules/videostab/src/wobble_suppression.cpp b/modules/videostab/src/wobble_suppression.cpp index 450ec0c441..a09437fa39 100644 --- a/modules/videostab/src/wobble_suppression.cpp +++ b/modules/videostab/src/wobble_suppression.cpp @@ -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); } diff --git a/samples/cpp/videostab.cpp b/samples/cpp/videostab.cpp index f401dccc55..dccd413066 100644 --- a/samples/cpp/videostab.cpp +++ b/samples/cpp/videostab.cpp @@ -86,6 +86,18 @@ void printHelp() " --stdev=(|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=(|auto)\n" + " Trimming ratio used in linear programming based method.\n" + " --lp-w1=(|1)\n" + " 1st derivative weight. The default is 1.\n" + " --lp-w2=(|10)\n" + " 2nd derivative weight. The default is 10.\n" + " --lp-w3=(|100)\n" + " 3rd derivative weight. The default is 100.\n" + " --lp-w4=(|100)\n" + " Non-translation motion components weight. The default is 100.\n\n" " --deblur=(yes|no)\n" " Do deblurring.\n" " --deblur-sens=\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(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(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;