added gain compensation into opencv_stitching

modules/stitching/blenders.cpp

@@ -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);

modules/stitching/blenders.hpp

@@ -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);

modules/stitching/exposure_compensate.cpp

@@ -45,3 +45,88 @@
 
 using namespace std;
 using namespace cv;
+
+
+Ptr<ExposureCompensator> ExposureCompensator::createDefault(int type)
+{
+    if (type == NO)
+        return new NoExposureCompensator();
+    if (type == OVERLAP)
+        return new OverlapExposureCompensator();
+    CV_Error(CV_StsBadArg, "unsupported exposure compensation method");
+    return NULL;
+}
+
+
+void OverlapExposureCompensator::feed(const vector<Point> &corners, const vector<Mat> &images, 
+                                      const vector<Mat> &masks)
+{
+    const int num_images = static_cast<int>(images.size());
+    Mat_<double> N(num_images, num_images); N.setTo(0);
+    Mat_<double> I(num_images, num_images); I.setTo(0);
+
+    Rect dst_roi = resultRoi(corners, images);
+    Mat subimg1, subimg2;
+    Mat_<uchar> submask1, submask2, overlap;
+
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = i; j < num_images; ++j)
+        {
+            Rect roi;
+            if (overlapRoi(corners[i], corners[j], images[i].size(), images[j].size(), roi))
+            {
+                subimg1 = images[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));
+                subimg2 = images[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));
+
+                submask1 = masks[i](Rect(roi.tl() - corners[i], roi.br() - corners[i]));
+                submask2 = masks[j](Rect(roi.tl() - corners[j], roi.br() - corners[j]));
+                overlap = submask1 & submask2;
+
+                N(i, j) = N(j, i) = countNonZero(overlap);
+
+                double Isum1 = 0, Isum2 = 0;
+                for (int y = 0; y < roi.height; ++y)
+                {
+                    const Point3_<uchar>* r1 = subimg1.ptr<Point3_<uchar> >(y);
+                    const Point3_<uchar>* r2 = subimg2.ptr<Point3_<uchar> >(y);
+                    for (int x = 0; x < roi.width; ++x)
+                    {
+                        if (overlap(y, x))
+                        {
+                            Isum1 += sqrt(static_cast<double>(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z)));
+                            Isum2 += sqrt(static_cast<double>(sqr(r2[x].x) + sqr(r2[x].y) + sqr(r2[x].z)));
+                        }
+                    }
+                }
+                I(i, j) = Isum1 / N(i, j);
+                I(j, i) = Isum2 / N(i, j);
+            }
+        }
+    }
+
+    double alpha = 0.01;
+    double beta = 100;
+
+    Mat_<double> A(num_images, num_images); A.setTo(0);
+    Mat_<double> b(num_images, 1); b.setTo(0);
+    for (int i = 0; i < num_images; ++i)
+    {
+        for (int j = 0; j < num_images; ++j)
+        {
+            b(i, 0) += beta * N(i, j);
+            A(i, i) += beta * N(i, j);
+            if (j == i) continue;
+            A(i, i) += 2 * alpha * I(i, j) * I(i, j) * N(i, j);
+            A(i, j) -= 2 * alpha * I(i, j) * I(j, i) * N(i, j);
+        }
+    }
+
+    solve(A, b, gains_, DECOMP_SVD);
+}
+
+
+void OverlapExposureCompensator::apply(int index, Point /*corner*/, Mat &image, const Mat &/*mask*/)
+{
+    image *= gains_(index, 0);
+}

modules/stitching/exposure_compensate.hpp

@@ -48,16 +48,33 @@
 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_;
 };
 
 

modules/stitching/main.cpp

@@ -39,6 +39,13 @@
 // the use of this software, even if advised of the possibility of such damage.
 //
 //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"
 #include "util.hpp"
 #include "warpers.hpp"
@@ -63,6 +70,7 @@ void printUsage()
         << "\t[--conf_thresh <float>]\n"
         << "\t[--wavecorrect (no|yes)]\n"
         << "\t[--warp (plane|cylindrical|spherical)]\n" 
+        << "\t[--exposcomp (no|overlap)]\n"
         << "\t[--seam (no|voronoi|graphcut)]\n" 
         << "\t[--blend (no|feather|multiband)]\n"
         << "\t[--numbands <int>]\n"
@@ -84,6 +92,7 @@ int ba_space = BundleAdjuster::FOCAL_RAY_SPACE;
 float conf_thresh = 1.f;
 bool wave_correct = true;
 int warp_type = Warper::SPHERICAL;
+int expos_comp_type = ExposureCompensator::OVERLAP;
 bool user_match_conf = false;
 float match_conf = 0.6f;
 int seam_find_type = SeamFinder::GRAPH_CUT;
@@ -186,6 +195,19 @@ int parseCmdArgs(int argc, char** argv)
             }
             i++;
         }
+        else if (string(argv[i]) == "--exposcomp")
+        {
+            if (string(argv[i + 1]) == "no")
+                expos_comp_type = ExposureCompensator::NO;
+            else if (string(argv[i + 1]) == "overlap")
+                expos_comp_type = ExposureCompensator::OVERLAP;
+            else
+            {
+                cout << "Bad exposure compensation method\n";
+                return -1;
+            }
+            i++;
+        }        
         else if (string(argv[i]) == "--seam")
         {
             if (string(argv[i + 1]) == "no")
@@ -372,7 +394,7 @@ int main(int argc, char* argv[])
         focals.push_back(cameras[i].focal);
     }
     nth_element(focals.begin(), focals.end(), focals.begin() + focals.size() / 2);
-    float camera_focal = static_cast<float>(focals[focals.size() / 2]);
+    float warped_image_scale = static_cast<float>(focals[focals.size() / 2]);
 
     LOGLN("Warping images (auxiliary)... ");
     t = getTickCount();
@@ -391,7 +413,7 @@ int main(int argc, char* argv[])
     }
 
     // Warp images and their masks
-    Ptr<Warper> warper = Warper::createByCameraFocal(static_cast<float>(camera_focal * seam_work_aspect), 
+    Ptr<Warper> warper = Warper::createByCameraFocal(static_cast<float>(warped_image_scale * seam_work_aspect), 
                                                      warp_type);
     for (int i = 0; i < num_images; ++i)
     {
@@ -410,8 +432,8 @@ int main(int argc, char* argv[])
 
     LOGLN("Exposure compensation (feed)...");
     t = getTickCount();
-    Ptr<ExposureCompensator> compensator = new NoExposureCompensator();
-    compensator->feed(images_warped, masks_warped);
+    Ptr<ExposureCompensator> compensator = ExposureCompensator::createDefault(expos_comp_type);
+    compensator->feed(corners, images_warped, masks_warped);
     LOGLN("Exposure compensation (feed), time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
 
     LOGLN("Finding seams...");
@@ -446,10 +468,27 @@ int main(int argc, char* argv[])
             if (compose_megapix > 0)
                 compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img.size().area()));
             is_compose_scale_set = true;
+
+            // Compute relative scales
             compose_seam_aspect = compose_scale / seam_scale;
             compose_work_aspect = compose_scale / work_scale;
-            camera_focal *= static_cast<float>(compose_work_aspect);
-            warper = Warper::createByCameraFocal(camera_focal, warp_type);
+
+            // Update warped image scale
+            warped_image_scale *= static_cast<float>(compose_work_aspect);
+            warper = Warper::createByCameraFocal(warped_image_scale, warp_type);
+
+            // Update corners and sizes
+            Rect dst_roi = resultRoi(corners, sizes);
+            for (int i = 0; i < num_images; ++i)
+            {
+                // Update camera focal
+                cameras[i].focal *= compose_work_aspect;
+
+                // Update corner and size
+                corners[i] = dst_roi.tl() + (corners[i] - dst_roi.tl()) * compose_seam_aspect;
+                sizes[i] = Size(static_cast<int>((sizes[i].width + 1) * compose_seam_aspect), 
+                                static_cast<int>((sizes[i].height + 1) * compose_seam_aspect));
+            }
         }
         if (abs(compose_scale - 1) > 1e-1)
             resize(full_img, img, Size(), compose_scale, compose_scale);
@@ -458,21 +497,18 @@ int main(int argc, char* argv[])
         full_img.release();
         Size img_size = img.size();
 
-        // Update cameras paramters
-        cameras[img_idx].focal *= compose_work_aspect;
-
         // Warp the current image
         warper->warp(img, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R, 
                      img_warped);
 
-        // Warp current image mask
+        // Warp the current image mask
         mask.create(img_size, CV_8U);
         mask.setTo(Scalar::all(255));    
         warper->warp(mask, static_cast<float>(cameras[img_idx].focal), cameras[img_idx].R, mask_warped,
                      INTER_NEAREST, BORDER_CONSTANT);
 
         // Compensate exposure
-        compensator->apply(img_idx, img_warped, mask_warped);
+        compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped);
 
         img_warped.convertTo(img_warped_s, CV_16S);
         img_warped.release();
@@ -486,21 +522,11 @@ int main(int argc, char* argv[])
         if (static_cast<Blender*>(blender) == 0)
         {
             blender = Blender::createDefault(blend_type);
-
             if (blend_type == Blender::MULTI_BAND)
             {
                 MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(static_cast<Blender*>(blender));
                 mb->setNumBands(numbands);
             }
-
-            // Determine the final image size
-            Rect dst_roi = resultRoi(corners, sizes);
-            for (int i = 0; i < num_images; ++i)
-            {
-                corners[i] = dst_roi.tl() + (corners[i] - dst_roi.tl()) * compose_seam_aspect;
-                sizes[i] = Size(static_cast<int>((sizes[i].width + 1) * compose_seam_aspect), 
-                                static_cast<int>((sizes[i].height + 1) * compose_seam_aspect));
-            }
             blender->prepare(corners, sizes);
         }
 

modules/stitching/seam_finders.cpp

@@ -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);
 

modules/stitching/util.cpp

@@ -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;
+}
+

modules/stitching/util.hpp

@@ -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__

modules/stitching/util_inl.hpp

@@ -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__