added gain compensation into opencv_stitching
authorAlexey Spizhevoy <no@email>
Wed, 25 May 2011 09:09:41 +0000 (09:09 +0000)
committerAlexey Spizhevoy <no@email>
Wed, 25 May 2011 09:09:41 +0000 (09:09 +0000)
modules/stitching/blenders.cpp
modules/stitching/blenders.hpp
modules/stitching/exposure_compensate.cpp
modules/stitching/exposure_compensate.hpp
modules/stitching/main.cpp
modules/stitching/seam_finders.cpp
modules/stitching/util.cpp
modules/stitching/util.hpp
modules/stitching/util_inl.hpp

index cce1381..916b5da 100644 (file)
@@ -282,24 +282,6 @@ void MultiBandBlender::blend(Mat &dst, Mat &dst_mask)
 //////////////////////////////////////////////////////////////////////////////
 // Auxiliary functions
 
-Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)
-{
-    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
-    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());
-
-    CV_Assert(sizes.size() == corners.size());
-    for (size_t i = 0; i < corners.size(); ++i)
-    {
-        tl.x = min(tl.x, corners[i].x);
-        tl.y = min(tl.y, corners[i].y);
-        br.x = max(br.x, corners[i].x + sizes[i].width);
-        br.y = max(br.y, corners[i].y + sizes[i].height);
-    }
-
-    return Rect(tl, br);
-}
-
-
 void normalize(const Mat& weight, Mat& src)
 {
     CV_Assert(weight.type() == CV_32F);
index ffd29e4..f65b8dd 100644 (file)
@@ -102,8 +102,6 @@ private:
 //////////////////////////////////////////////////////////////////////////////
 // Auxiliary functions
 
-cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Size> &sizes);
-
 void normalize(const cv::Mat& weight, cv::Mat& src);
 
 void createWeightMap(const cv::Mat& mask, float sharpness, cv::Mat& weight);
index c449159..db2cbc2 100644 (file)
 \r
 using namespace std;\r
 using namespace cv;\r
+\r
+\r
+Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)\r
+{\r
+    if (type == NO)
+        return new NoExposureCompensator();
+    if (type == OVERLAP)
+        return new OverlapExposureCompensator();
+    CV_Error(CV_StsBadArg, "unsupported exposure compensation method");
+    return NULL;\r
+}\r
+\r
+\r
+void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector<Mat> &images, \r
+                                      const vector<Mat> &masks)\r
+{\r
+    const int num_images = static_cast<int>(images.size());\r
+    Mat_<double> N(num_images, num_images); N.setTo(0);\r
+    Mat_<double> I(num_images, num_images); I.setTo(0);\r
+\r
+    Rect dst_roi = resultRoi(corners, images);\r
+    Mat subimg1, subimg2;\r
+    Mat_<uchar> submask1, submask2, overlap;\r
+\r
+    for (int i = 0; i < num_images; ++i)\r
+    {\r
+        for (int j = i; j < num_images; ++j)\r
+        {\r
+            Rect roi;\r
+            if (overlapRoi(corners[i], corners[j], images[i].size(), images[j].size(), roi))\r
+            {\r
+                subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));\r
+                subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));\r
+\r
+                submask1 = masks[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));\r
+                submask2 = masks[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));\r
+                overlap = submask1 & submask2;\r
+\r
+                N(i, j) = N(j, i) = countNonZero(overlap);\r
+\r
+                double Isum1 = 0, Isum2 = 0;\r
+                for (int y = 0; y < roi.height; ++y)\r
+                {\r
+                    const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(y);\r
+                    const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);\r
+                    for (int x = 0; x < roi.width; ++x)\r
+                    {\r
+                        if (overlap(y, x))\r
+                        {\r
+                            Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));\r
+                            Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));\r
+                        }\r
+                    }\r
+                }\r
+                I(i, j) = Isum1 / N(i, j);\r
+                I(j, i) = Isum2 / N(i, j);\r
+            }\r
+        }\r
+    }\r
+\r
+    double alpha = 0.01;\r
+    double beta = 100;\r
+\r
+    Mat_<double> A(num_images, num_images); A.setTo(0);\r
+    Mat_<double> b(num_images, 1); b.setTo(0);\r
+    for (int i = 0; i < num_images; ++i)\r
+    {\r
+        for (int j = 0; j < num_images; ++j)\r
+        {\r
+            b(i, 0) += beta * N(i, j);\r
+            A(i, i) += beta * N(i, j);\r
+            if (j == i) continue;\r
+            A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);\r
+            A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);\r
+        }\r
+    }\r
+\r
+    solve(A, b, gains_, DECOMP_SVD);\r
+}\r
+\r
+\r
+void OverlapExposureCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)\r
+{\r
+    image *= gains_(index, 0);\r
+}\r
index 9793854..5aec7c7 100644 (file)
 class ExposureCompensator
 {
 public:
-    virtual void feed(const std::vector<cv::Mat> &images, const std::vector<cv::Mat> &masks) = 0;
-    virtual void apply(int index, cv::Mat &image, cv::Mat &mask) = 0;
+    enum { NO, OVERLAP };
+    static cv::Ptr<ExposureCompensator> createDefault(int type);
+
+    virtual void feed(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images, 
+                      const std::vector<cv::Mat> &masks) = 0;
+    virtual void apply(int index, cv::Point corner, cv::Mat &image, const cv::Mat &mask) = 0;
 };
 
 
 class NoExposureCompensator : public ExposureCompensator
 {
 public:
-    virtual void feed(const std::vector<cv::Mat> &/*images*/, const std::vector<cv::Mat> &/*masks*/) {};
-    virtual void apply(int /*index*/, cv::Mat &/*image*/, cv::Mat &/*mask*/) {};
+    void feed(const std::vector<cv::Point> &/*corners*/, const std::vector<cv::Mat> &/*images*/, 
+              const std::vector<cv::Mat> &/*masks*/) {};
+    void apply(int /*index*/, cv::Point /*corner*/, cv::Mat &/*image*/, const cv::Mat &/*mask*/) {};
+};
+
+
+class OverlapExposureCompensator : public ExposureCompensator
+{
+public:
+    void feed(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images, 
+              const std::vector<cv::Mat> &masks);
+    void apply(int index, cv::Point corner, cv::Mat &image, const cv::Mat &mask);
+
+private:
+    cv::Mat_<double> gains_;
 };
 
 
index 9b47ca5..3856c8d 100644 (file)
 // the use of this software, even if advised of the possibility of such damage.\r
 //\r
 //M*/
+
+// We follow to methods described in these two papers:
+// - Heung-Yeung Shum and Richard Szeliski. 
+//   Construction of panoramic mosaics with global and local alignment. 2000.
+// - Matthew Brown and David G. Lowe. 
+//   Automatic Panoramic Image Stitching using Invariant Features. 2007.
+
 #include "precomp.hpp"\r
 #include "util.hpp"\r
 #include "warpers.hpp"\r
@@ -63,6 +70,7 @@ void printUsage()
         << "\t[--conf_thresh <float>]\n"\r
         << "\t[--wavecorrect (no|yes)]\n"\r
         << "\t[--warp (plane|cylindrical|spherical)]\n" \r
+        << "\t[--exposcomp (no|overlap)]\n"\r
         << "\t[--seam (no|voronoi|graphcut)]\n" \r
         << "\t[--blend (no|feather|multiband)]\n"\r
         << "\t[--numbands <int>]\n"\r
@@ -84,6 +92,7 @@ int ba_space = BundleAdjuster::FOCAL_RAY_SPACE;
 float conf_thresh = 1.f;\r
 bool wave_correct = true;\r
 int warp_type = Warper::SPHERICAL;\r
+int expos_comp_type = ExposureCompensator::OVERLAP;\r
 bool user_match_conf = false;\r
 float match_conf = 0.6f;\r
 int seam_find_type = SeamFinder::GRAPH_CUT;\r
@@ -186,6 +195,19 @@ int parseCmdArgs(int argc, char** argv)
             }\r
             i++;\r
         }\r
+        else if (string(argv[i]) == "--exposcomp")\r
+        {\r
+            if (string(argv[i + 1]) == "no")\r
+                expos_comp_type = ExposureCompensator::NO;\r
+            else if (string(argv[i + 1]) == "overlap")\r
+                expos_comp_type = ExposureCompensator::OVERLAP;\r
+            else\r
+            {\r
+                cout << "Bad exposure compensation method\n";\r
+                return -1;\r
+            }\r
+            i++;\r
+        }        \r
         else if (string(argv[i]) == "--seam")\r
         {\r
             if (string(argv[i + 1]) == "no")\r
@@ -372,7 +394,7 @@ int main(int argc, char* argv[])
         focals.push_back(cameras[i].focal);\r
     }\r
     nth_element(focals.begin(), focals.end(), focals.begin() + focals.size() / 2);\r
-    float camera_focal = static_cast<float>(focals[focals.size() / 2]);\r
+    float warped_image_scale = static_cast<float>(focals[focals.size() / 2]);\r
 \r
     LOGLN("Warping images (auxiliary)... ");\r
     t = getTickCount();\r
@@ -391,7 +413,7 @@ int main(int argc, char* argv[])
     }\r
 \r
     // Warp images and their masks\r
-    Ptr<Warper> warper = Warper::createByCameraFocal(static_cast<float>(camera_focal * seam_work_aspect), \r
+    Ptr<Warper> warper = Warper::createByCameraFocal(static_cast<float>(warped_image_scale * seam_work_aspect), \r
                                                      warp_type);\r
     for (int i = 0; i < num_images; ++i)\r
     {\r
@@ -410,8 +432,8 @@ int main(int argc, char* argv[])
 \r
     LOGLN("Exposure compensation (feed)...");\r
     t = getTickCount();\r
-    Ptr<ExposureCompensator> compensator = new NoExposureCompensator();\r
-    compensator->feed(images_warped, masks_warped);\r
+    Ptr<ExposureCompensator> compensator = ExposureCompensator::createDefault(expos_comp_type);\r
+    compensator->feed(corners, images_warped, masks_warped);\r
     LOGLN("Exposure compensation (feed), time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");\r
 \r
     LOGLN("Finding seams...");\r
@@ -446,10 +468,27 @@ int main(int argc, char* argv[])
             if (compose_megapix > 0)\r
                 compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img.size().area()));\r
             is_compose_scale_set = true;\r
+\r
+            // Compute relative scales\r
             compose_seam_aspect = compose_scale / seam_scale;\r
             compose_work_aspect = compose_scale / work_scale;\r
-            camera_focal *= static_cast<float>(compose_work_aspect);\r
-            warper = Warper::createByCameraFocal(camera_focal, warp_type);\r
+\r
+            // Update warped image scale\r
+            warped_image_scale *= static_cast<float>(compose_work_aspect);\r
+            warper = Warper::createByCameraFocal(warped_image_scale, warp_type);\r
+\r
+            // Update corners and sizes\r
+            Rect dst_roi = resultRoi(corners, sizes);\r
+            for (int i = 0; i < num_images; ++i)\r
+            {\r
+                // Update camera focal\r
+                cameras[i].focal *= compose_work_aspect;\r
+\r
+                // Update corner and size\r
+                corners[i] = dst_roi.tl() + (corners[i] - dst_roi.tl()) * compose_seam_aspect;\r
+                sizes[i] = Size(static_cast<int>((sizes[i].width + 1) * compose_seam_aspect), \r
+                                static_cast<int>((sizes[i].height + 1) * compose_seam_aspect));\r
+            }\r
         }\r
         if (abs(compose_scale - 1) > 1e-1)\r
             resize(full_img, img, Size(), compose_scale, compose_scale);\r
@@ -458,21 +497,18 @@ int main(int argc, char* argv[])
         full_img.release();\r
         Size img_size = img.size();\r
 \r
-        // Update cameras paramters\r
-        cameras[img_idx].focal *= compose_work_aspect;\r
-\r
         // Warp the current image\r
         warper->warp(img, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R, \r
                      img_warped);\r
 \r
-        // Warp current image mask\r
+        // Warp the current image mask\r
         mask.create(img_size, CV_8U);\r
         mask.setTo(Scalar::all(255));    \r
         warper->warp(mask, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R, mask_warped,\r
                      INTER_NEAREST, BORDER_CONSTANT);\r
 \r
         // Compensate exposure\r
-        compensator->apply(img_idx, img_warped, mask_warped);\r
+        compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped);\r
 \r
         img_warped.convertTo(img_warped_s, CV_16S);\r
         img_warped.release();\r
@@ -486,21 +522,11 @@ int main(int argc, char* argv[])
         if (static_cast<Blender*>(blender) == 0)\r
         {\r
             blender = Blender::createDefault(blend_type);\r
-\r
             if (blend_type == Blender::MULTI_BAND)\r
             {\r
                 MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(static_cast<Blender*>(blender));\r
                 mb->setNumBands(numbands);\r
             }\r
-\r
-            // Determine the final image size\r
-            Rect dst_roi = resultRoi(corners, sizes);\r
-            for (int i = 0; i < num_images; ++i)\r
-            {\r
-                corners[i] = dst_roi.tl() + (corners[i] - dst_roi.tl()) * compose_seam_aspect;\r
-                sizes[i] = Size(static_cast<int>((sizes[i].width + 1) * compose_seam_aspect), \r
-                                static_cast<int>((sizes[i].height + 1) * compose_seam_aspect));\r
-            }\r
             blender->prepare(corners, sizes);\r
         }\r
 \r
index eca665b..89f7394 100644 (file)
@@ -64,20 +64,13 @@ void PairwiseSeamFinder::find(const vector<Mat> &src, const vector<Point> &corne
 {
     if (src.size() == 0)
         return;
-
     for (size_t i = 0; i < src.size() - 1; ++i)
     {
         for (size_t j = i + 1; j < src.size(); ++j)
         {
-            int x_min = max(corners[i].x, corners[j].x);
-            int x_max = min(corners[i].x + src[i].cols - 1, corners[j].x + src[j].cols - 1);
-            int y_min = max(corners[i].y, corners[j].y);
-            int y_max = min(corners[i].y + src[i].rows - 1, corners[j].y + src[j].rows - 1);
-
-            if (x_max >= x_min && y_max >= y_min)
-                findInPair(src[i], src[j], corners[i], corners[j],
-                           Rect(x_min, y_min, x_max - x_min + 1, y_max - y_min + 1),
-                           masks[i], masks[j]);
+            Rect roi;
+            if (overlapRoi(corners[i], corners[j], src[i].size(), src[j].size(), roi))
+                findInPair(src[i], src[j], corners[i], corners[j], roi, masks[i], masks[j]);
         }
     }
 }
@@ -87,7 +80,6 @@ void VoronoiSeamFinder::findInPair(const Mat &img1, const Mat &img2, Point tl1,
                                    Rect roi, Mat &mask1, Mat &mask2)
 {
     const int gap = 10;
-
     Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
     Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U);
 
index 88162dd..1bcb7f8 100644 (file)
@@ -44,7 +44,8 @@
 using namespace std;
 using namespace cv;
 
-void DjSets::create(int n) {
+void DjSets::create(int n) 
+{
     rank_.assign(n, 0);
     size.assign(n, 1);
     parent.resize(n);
@@ -53,12 +54,14 @@ void DjSets::create(int n) {
 }
 
 
-int DjSets::find(int elem) {
+int DjSets::find(int elem) 
+{
     int set = elem;
     while (set != parent[set])
         set = parent[set];
     int next;
-    while (elem != parent[elem]) {
+    while (elem != parent[elem]) 
+    {
         next = parent[elem];
         parent[elem] = set;
         elem = next;
@@ -67,13 +70,16 @@ int DjSets::find(int elem) {
 }
 
 
-int DjSets::merge(int set1, int set2) {
-    if (rank_[set1] < rank_[set2]) {
+int DjSets::merge(int set1, int set2) 
+{
+    if (rank_[set1] < rank_[set2]) 
+    {
         parent[set1] = set2;
         size[set2] += size[set1];
         return set2;
     }
-    if (rank_[set2] < rank_[set1]) {
+    if (rank_[set2] < rank_[set1]) 
+    {
         parent[set2] = set1;
         size[set1] += size[set2];
         return set1;
@@ -89,3 +95,56 @@ void Graph::addEdge(int from, int to, float weight)
 {
     edges_[from].push_back(GraphEdge(from, to, weight));
 }
+
+
+bool overlapRoi(Point tl1, Point tl2, Size sz1, Size sz2, Rect &roi)
+{
+    int x_tl = max(tl1.x, tl2.x);
+    int y_tl = max(tl1.y, tl2.y);
+    int x_br = min(tl1.x + sz1.width, tl2.x + sz2.width);
+    int y_br = min(tl1.y + sz1.height, tl2.y + sz2.height);
+    if (x_tl < x_br && y_tl < y_br)
+    {
+        roi = Rect(x_tl, y_tl, x_br - x_tl, y_br - y_tl);
+        return true;
+    }
+    return false;
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Mat> &images)
+{
+    vector<Size> sizes(images.size());
+    for (size_t i = 0; i < images.size(); ++i)
+        sizes[i] = images[i].size();
+    return resultRoi(corners, sizes);
+}
+
+
+Rect resultRoi(const vector<Point> &corners, const vector<Size> &sizes)
+{
+    CV_Assert(sizes.size() == corners.size());
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    Point br(numeric_limits<int>::min(), numeric_limits<int>::min());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+        br.x = max(br.x, corners[i].x + sizes[i].width);
+        br.y = max(br.y, corners[i].y + sizes[i].height);
+    }
+    return Rect(tl, br);
+}
+
+
+Point resultTl(const vector<Point> &corners)
+{
+    Point tl(numeric_limits<int>::max(), numeric_limits<int>::max());
+    for (size_t i = 0; i < corners.size(); ++i)
+    {
+        tl.x = min(tl.x, corners[i].x);
+        tl.y = min(tl.y, corners[i].y);
+    }
+    return tl;
+}
+
index 60ac074..564fc25 100644 (file)
@@ -90,23 +90,26 @@ class Graph
 {
 public:
     Graph(int num_vertices = 0) { create(num_vertices); }
-
     void create(int num_vertices) { edges_.assign(num_vertices, std::list<GraphEdge>()); }
-
     int numVertices() const { return static_cast<int>(edges_.size()); }
-
     void addEdge(int from, int to, float weight);
-
-    template <typename B>
-    B forEach(B body) const;
-
-    template <typename B> 
-    B walkBreadthFirst(int from, B body) const;
+    template <typename B> B forEach(B body) const;
+    template <typename B> B walkBreadthFirst(int from, B body) const;
     
 private:
     std::vector< std::list<GraphEdge> > edges_;
 };
 
+
+//////////////////////////////////////////////////////////////////////////////
+// Auxiliary functions
+
+bool overlapRoi(cv::Point tl1, cv::Point tl2, cv::Size sz1, cv::Size sz2, cv::Rect &roi);
+cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Mat> &images);
+cv::Rect resultRoi(const std::vector<cv::Point> &corners, const std::vector<cv::Size> &sizes);
+cv::Point resultTl(const std::vector<cv::Point> &corners);
+
+
 #include "util_inl.hpp"
 
 #endif // __OPENCV_STITCHING_UTIL_HPP__
index bdd1176..12f8396 100644 (file)
@@ -105,11 +105,8 @@ double normL2sq(const cv::Mat &r)
 }
 
 
-template <typename T>
-static inline
-T sqr(T x)
-{
-    return x * x;
-}
+static inline int sqr(int x) { return x * x; }
+static inline float sqr(float x) { return x * x; }
+static inline double sqr(double x) { return x * x; }
 
 #endif // __OPENCV_STITCHING_UTIL_INL_HPP__