diff --git a/modules/xphoto/CMakeLists.txt b/modules/xphoto/CMakeLists.txt index 3e550a3c892..3a2a170988c 100644 --- a/modules/xphoto/CMakeLists.txt +++ b/modules/xphoto/CMakeLists.txt @@ -1,2 +1,2 @@ set(the_description "Addon to basic photo module") -ocv_define_module(xphoto opencv_core opencv_imgproc WRAP python) +ocv_define_module(xphoto opencv_core opencv_imgproc opencv_optflow opencv_ximgproc opencv_calib3d opencv_highgui WRAP python) diff --git a/modules/xphoto/README.md b/modules/xphoto/README.md index fc4aa8857b9..217ee269332 100644 --- a/modules/xphoto/README.md +++ b/modules/xphoto/README.md @@ -4,4 +4,5 @@ Additional photo processing algorithms 1. Color balance 2. Denoising 3. Inpainting +4. Obstruction-free photography (occlusion and reflection) diff --git a/modules/xphoto/doc/xphoto.bib b/modules/xphoto/doc/xphoto.bib index b57471ceb64..9d6ce705c8e 100644 --- a/modules/xphoto/doc/xphoto.bib +++ b/modules/xphoto/doc/xphoto.bib @@ -14,3 +14,12 @@ @inproceedings{Cheng2015 pages={1000--1008}, year={2015} } + +@article{Xue2015ObstructionFree, + author = {Tianfan Xue and Michael Rubinstein and Ce Liu and William T. Freeman}, + title = {A Computational Approach for Obstruction-Free Photography}, + journal = {ACM Transactions on Graphics (Proc. SIGGRAPH)}, + year = {2015}, + volume = {34}, + number = {4}, +} diff --git a/modules/xphoto/include/opencv2/xphoto.hpp b/modules/xphoto/include/opencv2/xphoto.hpp index 73b1e0bb3c4..70b35a1fd9c 100644 --- a/modules/xphoto/include/opencv2/xphoto.hpp +++ b/modules/xphoto/include/opencv2/xphoto.hpp @@ -50,4 +50,5 @@ #include "xphoto/white_balance.hpp" #include "xphoto/dct_image_denoising.hpp" #include "xphoto/bm3d_image_denoising.hpp" +#include "xphoto/obstruction_free.hpp" #endif diff --git a/modules/xphoto/include/opencv2/xphoto/obstruction_free.hpp b/modules/xphoto/include/opencv2/xphoto/obstruction_free.hpp new file mode 100644 index 00000000000..ff510a2b070 --- /dev/null +++ b/modules/xphoto/include/opencv2/xphoto/obstruction_free.hpp @@ -0,0 +1,157 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + + +#ifndef __OPENCV_OBSTRUCTION_FREE_HPP__ +#define __OPENCV_OBSTRUCTION_FREE_HPP__ + +/** @file +@date June 18, 2017 +@author Binbin Xu +*/ + +#include "opencv2/core.hpp" +#include "opencv2/imgproc.hpp" +#include "opencv2/ximgproc/sparse_match_interpolator.hpp" +#include "opencv2/optflow.hpp" +#include "opencv2/calib3d.hpp" + +#include "opencv2/highgui.hpp" + +namespace cv +{ +namespace xphoto +{ + +//! @addtogroup xphoto +//! @{ + + + /** @brief The class implements a general obstruction free approach that can remove occlusions and reflections from input image sequences without manual masks. + + See the original paper @cite Xue2015ObstructionFree for more details. + + + */ + class CV_EXPORTS obstructionFree + { + public: + /*! + * @brief Constructors + */ + obstructionFree(); + + /*! + * @brief Constructors + * @param srcImgs input image sequences + */ + obstructionFree(const std::vector &srcImgs); + + /*! + * @brief core function to remove occlusions + * @param srcImgs source image sequences, involving translation motions. + * @param dst Obstruction-removed destination image, corresponding to the reference image, with the same size and type. In general, the reference image is the center frame of the input image. + * @param foreground estimated reflection or opaque obstruction layer + * @param mask estimated occlusion areas (CV_8UC1), where zero pixels indicate area to be estimated to be occlusions. + * @param obstructionType: reflection (0) or opaque obstruction (1) + */ + void removeOcc(const std::vector &srcImgs, Mat &dst, Mat& foreground, Mat &mask, const int obstructionType); + + private: + /*! + * @brief Parameters + */ + size_t frameNumber; //frame number of the input sequences + size_t referenceNumber; //target frame + int pyramidLevel; // Pyramid level + int coarseIterations; // iteration number for the coarsest level + int upperIterations; // iteration number for the upper level + int fixedPointIterations; // during each level of the pyramid + int sorIterations; // iterations of SOR + float omega; // relaxation factor in SOR + float lambda1; // weight for alpha map smoothness constraints + float lambda2; // weight for image smoothness constraints + float lambda3; // weight for independence between back/foreground component + float lambda4; // weight for gradient sparsity + + + std::vector backFlowFields; //estimated optical flow fields in the background layer + std::vector foreFlowFields; //estimated optical flow fields in the foreground layer + //std::vector warpedToReference; //warped images from input images through the estimated background flow fields + + //switch different methods. TODO + //int interpolationType; //interpolation type: 0(reflection) or 1(opaque occlusion) + //int edgeflowType; //edge flow type + + /** @brief private functions + */ + /** @brief Build pyramid by stacking all input image sequences + */ + std::vector buildPyramid( const std::vector & srcImgs); + + /** @brief Extract certain level of image sequences from the stacked image pyramid + */ + std::vector extractLevelImgs(const std::vector& pyramid, const int level); + + /** @brief Initialization: decompose the motion fields + */ + void motionInitDirect(const std::vector& video_input, std::vector& back_Flowfields, std::vector& fore_flowfields, std::vector& warpedToReference); + + /** @brief Initialization: decompose the image components in the case of reflection + */ + Mat imgInitDecomRef(const std::vector &warpedImgs); + + /** @brief Initialization: decompose the image components in the case of opaque reflection + */ + Mat imgInitDecomOpaq(const std::vector &warpedImgs, Mat& foreground, Mat& alphaMap); + + /** @brief Convert from sparse edge displacement to dense motion fields + */ + Mat sparseToDense(const Mat& im1, const Mat& im2, const Mat& im1_edges, const Mat& sparseFlow); + + /** @brief Visualize the optical flow flow with the window named figName + */ + void colorFlow(const Mat& flow, std::string figName); + + /** @brief Decompose motion between two images: target frame and source frame, using homography ransac + */ + void initMotionDecompose(const Mat& im1 , const Mat& im2 , Mat& back_denseFlow, Mat& fore_denseFlow, int back_ransacThre, int fore_ransacThre); + + /** @brief Warp im1 to output through optical flow: + */ + Mat imgWarpFlow(const Mat& im1, const Mat& flow); + + /** @brief Estimate homography matrix using edge flow fields + */ + Mat flowHomography(const Mat& edges, const Mat& flow, const int ransacThre); + + /** @brief Convert from index to matrix + */ + Mat indexToMask(const Mat& indexMat, const int rows, const int cols); + + /** @brief Calculate laplacian filters D_x^TD_x + D_y^TD_y + */ + Mat Laplac(const Mat& input); + + /*! + * @brief Calculate the weights in the alternative motion decomposition step + * @param inputSequence Input sequences + * @param background Estimated background component in the last iteration + * @param foreground Estimated foreground/obstruction component in the last iteration + * @param alphaMap Estimated alpha map for obstruction layer + * @return omega_1 The weight omega_1 = \phi(||I^t - W_O^t\hat{I}_O - W_O^t\hat{A} \circ W_B^t\hat{I}_B||^2)^{-1} + * @return omega_2 The weight omega_2 = \phi(||D_x\hat{I}_B||^2+||D_y\hat{I}_B||^2)^{-1} + * @return omega_3 The weight omega_3 = \phi(||D_x\hat{I}_O||^2+||D_y\hat{I}_O||^2)^{-1} + */ + void motDecomIrlsWeight(const std::vector& inputSequence, const Mat& background, const Mat& foreground, + Mat& alphaMap, std::vector& omega_1, std::vector& omega_2, std::vector& omega_3); + + }; +//! @} + + +} +} + +#endif // __OPENCV_OBSTRUCTION_FREE_HPP__ diff --git a/modules/xphoto/src/obstruction_free.cpp b/modules/xphoto/src/obstruction_free.cpp new file mode 100644 index 00000000000..14aab32419d --- /dev/null +++ b/modules/xphoto/src/obstruction_free.cpp @@ -0,0 +1,469 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html. + + + +#ifndef __OPENCV_OBSTRUCTION_FREE_CPP__ +#define __OPENCV_OBSTRUCTION_FREE_CPP__ + +#include +#include + +#include + +namespace cv +{ +namespace xphoto +{ + + +//constructor +obstructionFree::obstructionFree() +{ + //parameter + pyramidLevel = 3; // Pyramid level + coarseIterations = 4; // iteration number for the coarsest level + upperIterations = 1; // iteration number for the upper level + fixedPointIterations = 5; // during each level of the pyramid + sorIterations = 25; // iterations of SOR + omega = 1.6f; // relaxation factor in SOR + lambda1 = 1.0f; // weight for alpha map smoothness constraints + lambda2 = 0.1f; // weight for image smoothness constraints + lambda3 = 3000.0f; // weight for independence between back/foreground component + lambda4 = 0.5f; // weight for gradient sparsity +} + +obstructionFree::obstructionFree(const std::vector &srcImgs) +{ + //parameter + frameNumber = srcImgs.size(); + referenceNumber=(frameNumber-1)/2; //target frame + pyramidLevel = 3; // Pyramid level + coarseIterations = 4; // iteration number for the coarsest level + upperIterations = 1; // iteration number for the upper level + fixedPointIterations = 5; // during each level of the pyramid + sorIterations = 25; // iterations of SOR + omega = 1.6f; // relaxation factor in SOR + lambda1 = 1.0f; // weight for alpha map smoothness constraints + lambda2 = 0.1f; // weight for image smoothness constraints + lambda3 = 3000.0f; // weight for independence between back/foreground component + lambda4 = 0.5f; // weight for gradient sparsity +} + +//build pyramid by stacking all input image sequences +std::vector obstructionFree::buildPyramid(const std::vector & srcImgs) +{ + std::vector pyramid; + + for (size_t frame_i=0; frame_iframeNumber; frame_i++){ + Mat grey; + cvtColor(srcImgs[frame_i], grey, COLOR_RGB2GRAY); + pyramid.push_back(grey.clone()); + } + + Mat thisLevel, nextLevel; + size_t thisLevelId; + for (int level_i=1; level_ipyramidLevel; level_i++){ + for (size_t frame_i=0; frame_iframeNumber; frame_i++){ + thisLevelId=(level_i-1)*this->frameNumber+frame_i; + //nextLevelId=i*this.pyramidLevel+frame_i; + thisLevel=pyramid[thisLevelId]; + pyrDown(thisLevel, nextLevel); + pyramid.push_back(nextLevel.clone()); + } + } + return pyramid; +} + +//extract certain level of image sequences from the stacked image pyramid +std::vector obstructionFree::extractLevelImgs(const std::vector& pyramid, const int level){ + std::vector levelPyramid; + size_t imgId; + for (size_t frame_i=0; frame_iframeNumber; frame_i++){ + imgId=level*this->frameNumber+frame_i; + levelPyramid.push_back(pyramid[imgId].clone()); + } + return levelPyramid; +} + +Mat obstructionFree::indexToMask(const Mat& indexMat, const int rows, const int cols){ + Mat maskMat=Mat::zeros(rows, cols, CV_8UC1); + for (int i = 0; i < indexMat.cols; i++ ) { + for (int j = 0; j < indexMat.rows; j++) { + Vec2i mask_loca = indexMat.at(j, i); + if (mask_loca[0] !=0 && mask_loca[1] !=0) { + maskMat.at(Point(mask_loca)) = 255;} + }} + return maskMat; +} + +//estimate homography matrix using edge flow fields +Mat obstructionFree::flowHomography(const Mat& edges, const Mat& flow, const int ransacThre){ + Mat inlierMask, inlier_edges, inilier_edgeLocations; + std::vector edge_Locations; + + findNonZero(edges, edge_Locations); + + std::vector obj_edgeflow; + + for(size_t i = 0; i(src_y, src_x); + obj_edgeflow.push_back(Point2f(src_x + f.x, src_y + f.y)); + } + + Mat Homography = findHomography( edge_Locations, obj_edgeflow, RANSAC, ransacThre, inlierMask); + + Mat(edge_Locations).copyTo(inilier_edgeLocations,inlierMask); + + //convert index matrix to mask matrix + inlier_edges=indexToMask(inilier_edgeLocations, edges.rows, edges.cols); + + return inlier_edges; +} + +Mat obstructionFree::sparseToDense(const Mat& im1, const Mat& im2, const Mat& im1_edges, const Mat& sparseFlow){ + Mat denseFlow; + std::vector sparseFrom; + std::vector sparseTo; + + std::vector edge_Location; + findNonZero(im1_edges, edge_Location); + for(size_t i = 0; i(src_y, src_x); + sparseTo.push_back(Point2f(float(src_x + f.x), float(src_y + f.y))); + } + + Ptr epicInterpolation=ximgproc::createEdgeAwareInterpolator(); + epicInterpolation->interpolate(im1,sparseFrom,im2,sparseTo,denseFlow); + return denseFlow; +} + +//show motion fields using color +void obstructionFree::colorFlow(const Mat& flow, std::string figName="optical flow") +{ + //extraxt x and y channels + Mat xy[2]; //X,Y + split(flow, xy); + + //calculate angle and magnitude + Mat magnitude, angle; + cartToPolar(xy[0], xy[1], magnitude, angle, true); + + //translate magnitude to range [0;1] + double mag_max; + minMaxLoc(magnitude, 0, &mag_max); + magnitude.convertTo(magnitude, -1, 1.0 / mag_max); + + //build hsv image + Mat _hsv[3], hsv; + _hsv[0] = angle; + _hsv[1] = Mat::ones(angle.size(), CV_32F); + _hsv[2] = magnitude; + merge(_hsv, 3, hsv); + + //convert to BGR and show + Mat bgr;//CV_32FC3 matrix + cvtColor(hsv, bgr, COLOR_HSV2BGR); + imshow(figName, bgr); +} + +//motion decomposition between two images: target frame and source frame +void obstructionFree::initMotionDecompose(const Mat& im1, const Mat& im2, Mat& back_denseFlow, Mat& fore_denseFlow, int back_ransacThre=1, int fore_ransacThre=1){ + if (im1.channels()!= 1) + cvtColor(im1, im1, COLOR_RGB2GRAY); + if (im2.channels()!= 1) + cvtColor(im2, im2, COLOR_RGB2GRAY); + + Mat im1_edge, im2_edge; + Mat flow; + Mat edgeflow; //extracted edgeflow + + //Mat backH, mask_backH; + Mat back_edges, rest_edges, fore_edges; //edges aligned to the back layer using homography, remaining layer, foreground layers + Mat back_flow, rest_flow, fore_flow; + + + Canny(im1, im1_edge, 10, 100,3,true); + Canny(im2, im2_edge, 10, 100,3,true); + + ///////////////replace edgeflow + Ptr deepflow = optflow::createOptFlow_DeepFlow(); + deepflow->calc(im1, im2, flow); + //colorFlow(flow,"optical_flow"); + flow.copyTo(edgeflow, im1_edge); + //colorFlow(edgeflow,"edge_flow"); + +////////flow=>points using homography-ransac filtering, and then extract flow on the filtered edges + back_edges=flowHomography(im1_edge, edgeflow, back_ransacThre); + //imshow("back_edges", back_edges); + edgeflow.copyTo(back_flow,back_edges); + //colorFlow(back_flow, "back_flow"); + //////////rest edges and flows + rest_edges=im1_edge-back_edges; + //imshow("rest_edges", rest_edges); + rest_flow=edgeflow-back_flow; + //colorFlow(rest_flow, "rest_flow"); + + ////////////align resting flows to another homograghy + fore_edges=flowHomography(rest_edges, rest_flow, fore_ransacThre); + //imshow("fore_edges", fore_edges); + rest_flow.copyTo(fore_flow,fore_edges); + //colorFlow(fore_flow, "fore_flow"); + +///////////////////interpolation from sparse edgeFlow to denseFlow///////////////////// + back_denseFlow=sparseToDense(im1, im2, back_edges, back_flow); + fore_denseFlow=sparseToDense(im1, im2, fore_edges, fore_flow); + //colorFlow(back_denseFlow,"inter_back_denseflow"); + //colorFlow(fore_denseFlow,"inter_fore_denseflow"); +} + +//warping im1 to output through optical flow: +//flow=flow->cal(im1,im2), so warp im2 to back +Mat obstructionFree::imgWarpFlow(const Mat& im1, const Mat& flow){ + Mat flowmap_x(flow.size(), CV_32FC1); + Mat flowmap_y(flow.size(), CV_32FC1); + for (int j = 0; j < flowmap_x.rows; j++){ + for (int i = 0; i < flowmap_x.cols; ++i){ + Point2f f = flow.at(j, i); + flowmap_y.at(j, i) = float(j + f.y); + flowmap_x.at(j, i) = float(i + f.x); + }} + Mat warpedFrame; + remap(im1, warpedFrame, flowmap_x,flowmap_y ,INTER_CUBIC,BORDER_CONSTANT,255); + return warpedFrame; +} + +//Initialization: decompose the motion fields +void obstructionFree::motionInitDirect(const std::vector& video_input, std::vector& back_flowfields, std::vector& fore_flowfields, std::vector& warpedToReference){ + int back_ransacThre=1; + int fore_ransacThre=1; + + for (size_t frame_i=0; frame_i &warpedImgs){ + Mat background; + background=warpedImgs[referenceNumber]; + for (size_t frame_i=0; frame_i &warpedImgs, Mat& foreground, Mat& alphaMap){ + Mat background; + Mat sum=Mat::zeros(warpedImgs[referenceNumber].rows,warpedImgs[referenceNumber].cols,CV_32F); + Mat temp,background_temp; + for (size_t frame_i=0; frame_i(offset)=_input.at(offset+1) - _input.at(offset); + } + + for(int i=0;i(offset) -= temp.at(offset); + if(j>0) + _output.at(offset) += temp.at(offset-1); + } + + temp.release(); + + temp=Mat::zeros(s, CV_32FC1); + + // vertical filtering + for(int i=0;i(offset)=_input.at(offset+width)- _input.at(offset); + } + for(int i=0;i(offset) -= temp.at(offset); + if(i>0) + _output.at(offset) += temp.at(offset-width); + } + + output = _output.reshape(0, input.rows); + //output.convertTo(output, input.type()); + return output; + } + +//Calculate the weights in the alternative motion decomposition step +void obstructionFree::motDecomIrlsWeight(const std::vector& inputSequence, const Mat& background, const Mat& foreground, + Mat& alphaMap, std::vector& omega_1, std::vector& omega_2, std::vector& omega_3){ + + int width = background.cols; + int height = background.rows; + int npixels = width*height; + + + CV_Assert(omega_1.size() == 0 || omega_1.size() == inputSequence.size()*background.total()); + CV_Assert(omega_2.size() == 0 || omega_2.size() == background.total()); + CV_Assert(omega_3.size() == 0 || omega_3.size() == background.total()); + + float varepsilon = 0.000001f; + int deriv_ddepth = CV_32F; + + Mat backgd_dx; + Mat backgd_dy; + Mat obstruc_dx; + Mat obstruc_dy; + + + //compute gradients of current background and occlusion components + Sobel(background, backgd_dx, deriv_ddepth, 1, 0); + Sobel(background, backgd_dy, deriv_ddepth, 0, 1); + Sobel(foreground, obstruc_dx, deriv_ddepth, 1, 0); + Sobel(foreground, obstruc_dy, deriv_ddepth, 0, 1); + + //if weights have not been initialized + if (omega_1.size()==0 && omega_2.size() == 0 && omega_3.size() == 0){ + + //compute derivative denominators (weights) + for(size_t tt=0; ttbackFlowFields[tt]; + Mat obstruc_flow=this->foreFlowFields[tt]; + Mat temp = img-imgWarpFlow(foreground,obstruc_flow)-imgWarpFlow(alphaMap,obstruc_flow).mul(imgWarpFlow(background,back_flow)); + for(int i=0; i(i)*temp.at(i)+varepsilon)); + } + } + for(int i=0;i(i)*backgd_dx.at(i)+backgd_dy.at(i)*backgd_dy.at(i)+varepsilon)) ; + omega_3.push_back(1/sqrt(obstruc_dx.at(i)*obstruc_dx.at(i)+obstruc_dy.at(i)*obstruc_dy.at(i)+varepsilon)); + } + } + + //if weights have been calculated + else if (omega_1.size() == inputSequence.size()*background.total() && omega_2.size() == background.total() && + omega_3.size() == background.total()){ + + //compute derivative denominators (weights) + for(size_t tt=0; ttbackFlowFields[tt]; + Mat obstruc_flow=this->foreFlowFields[tt]; + Mat temp = img-imgWarpFlow(foreground,obstruc_flow)-imgWarpFlow(alphaMap,obstruc_flow).mul(imgWarpFlow(background,back_flow)); + for(int i=0; i(i + tt*npixels); + omega_1[offset] = 1/sqrt(temp.at(i)*temp.at(i)+varepsilon); + } + } + for(int i=0;i(i); + omega_2[offset] = 1/sqrt(backgd_dx.at(i)*backgd_dx.at(i)+backgd_dy.at(i)*backgd_dy.at(i)+varepsilon); + omega_3[offset] = 1/sqrt(obstruc_dx.at(i)*obstruc_dx.at(i)+obstruc_dy.at(i)*obstruc_dy.at(i)+varepsilon); + } + } +} + +//core function +void obstructionFree::removeOcc(const std::vector &srcImgs, Mat &dst, Mat& foreground, Mat &mask, const int obstructionType){ + //initialization + std::vector warpedToReference; + std::vector srcPyramid=buildPyramid(srcImgs); + std::vector coarseLevel = extractLevelImgs(srcPyramid,pyramidLevel-1); + motionInitDirect(coarseLevel, backFlowFields, foreFlowFields, warpedToReference); + + CV_Assert(obstructionType==0 || obstructionType==1); + if (obstructionType==0) + dst=imgInitDecomRef(warpedToReference); + else + dst=imgInitDecomOpaq(warpedToReference,foreground,mask); + + //alternative optimization:TODO + + +} + + + + +} +} +#endif // __OPENCV_OBSTRUCTION_FREE_CPP__