diff --git a/modules/stitching/blenders.cpp b/modules/stitching/blenders.cpp index cce1381b1a..916b5dae6e 100644 --- a/modules/stitching/blenders.cpp +++ b/modules/stitching/blenders.cpp @@ -282,24 +282,6 @@ void MultiBandBlender::blend(Mat &dst, Mat &dst_mask) ////////////////////////////////////////////////////////////////////////////// // Auxiliary functions -Rect resultRoi(const vector &corners, const vector &sizes) -{ - Point tl(numeric_limits::max(), numeric_limits::max()); - Point br(numeric_limits::min(), numeric_limits::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); diff --git a/modules/stitching/blenders.hpp b/modules/stitching/blenders.hpp index ffd29e4e27..f65b8dde94 100644 --- a/modules/stitching/blenders.hpp +++ b/modules/stitching/blenders.hpp @@ -102,8 +102,6 @@ private: ////////////////////////////////////////////////////////////////////////////// // Auxiliary functions -cv::Rect resultRoi(const std::vector &corners, const std::vector &sizes); - void normalize(const cv::Mat& weight, cv::Mat& src); void createWeightMap(const cv::Mat& mask, float sharpness, cv::Mat& weight); diff --git a/modules/stitching/exposure_compensate.cpp b/modules/stitching/exposure_compensate.cpp index c449159a1f..db2cbc20bf 100644 --- a/modules/stitching/exposure_compensate.cpp +++ b/modules/stitching/exposure_compensate.cpp @@ -45,3 +45,88 @@ using namespace std; using namespace cv; + + +Ptr 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 &corners, const vector &images, + const vector &masks) +{ + const int num_images = static_cast(images.size()); + Mat_ N(num_images, num_images); N.setTo(0); + Mat_ I(num_images, num_images); I.setTo(0); + + Rect dst_roi = resultRoi(corners, images); + Mat subimg1, subimg2; + Mat_ 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_* r1 = subimg1.ptr >(y); + const Point3_* r2 = subimg2.ptr >(y); + for (int x = 0; x < roi.width; ++x) + { + if (overlap(y, x)) + { + Isum1 += sqrt(static_cast(sqr(r1[x].x) + sqr(r1[x].y) + sqr(r1[x].z))); + Isum2 += sqrt(static_cast(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_ A(num_images, num_images); A.setTo(0); + Mat_ 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); +} diff --git a/modules/stitching/exposure_compensate.hpp b/modules/stitching/exposure_compensate.hpp index 9793854f41..5aec7c7d76 100644 --- a/modules/stitching/exposure_compensate.hpp +++ b/modules/stitching/exposure_compensate.hpp @@ -48,16 +48,33 @@ class ExposureCompensator { public: - virtual void feed(const std::vector &images, const std::vector &masks) = 0; - virtual void apply(int index, cv::Mat &image, cv::Mat &mask) = 0; + enum { NO, OVERLAP }; + static cv::Ptr createDefault(int type); + + virtual void feed(const std::vector &corners, const std::vector &images, + const std::vector &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 &/*images*/, const std::vector &/*masks*/) {}; - virtual void apply(int /*index*/, cv::Mat &/*image*/, cv::Mat &/*mask*/) {}; + void feed(const std::vector &/*corners*/, const std::vector &/*images*/, + const std::vector &/*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 &corners, const std::vector &images, + const std::vector &masks); + void apply(int index, cv::Point corner, cv::Mat &image, const cv::Mat &mask); + +private: + cv::Mat_ gains_; }; diff --git a/modules/stitching/main.cpp b/modules/stitching/main.cpp index 9b47ca567c..3856c8d5a9 100644 --- a/modules/stitching/main.cpp +++ b/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 ]\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 ]\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(focals[focals.size() / 2]); + float warped_image_scale = static_cast(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::createByCameraFocal(static_cast(camera_focal * seam_work_aspect), + Ptr warper = Warper::createByCameraFocal(static_cast(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 compensator = new NoExposureCompensator(); - compensator->feed(images_warped, masks_warped); + Ptr 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(compose_work_aspect); - warper = Warper::createByCameraFocal(camera_focal, warp_type); + + // Update warped image scale + warped_image_scale *= static_cast(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((sizes[i].width + 1) * compose_seam_aspect), + static_cast((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(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(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) == 0) { blender = Blender::createDefault(blend_type); - if (blend_type == Blender::MULTI_BAND) { MultiBandBlender* mb = dynamic_cast(static_cast(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((sizes[i].width + 1) * compose_seam_aspect), - static_cast((sizes[i].height + 1) * compose_seam_aspect)); - } blender->prepare(corners, sizes); } diff --git a/modules/stitching/seam_finders.cpp b/modules/stitching/seam_finders.cpp index eca665b168..89f73948ca 100644 --- a/modules/stitching/seam_finders.cpp +++ b/modules/stitching/seam_finders.cpp @@ -64,20 +64,13 @@ void PairwiseSeamFinder::find(const vector &src, const vector &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); diff --git a/modules/stitching/util.cpp b/modules/stitching/util.cpp index 88162ddc71..1bcb7f813f 100644 --- a/modules/stitching/util.cpp +++ b/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 &corners, const vector &images) +{ + vector 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 &corners, const vector &sizes) +{ + CV_Assert(sizes.size() == corners.size()); + Point tl(numeric_limits::max(), numeric_limits::max()); + Point br(numeric_limits::min(), numeric_limits::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 &corners) +{ + Point tl(numeric_limits::max(), numeric_limits::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; +} + diff --git a/modules/stitching/util.hpp b/modules/stitching/util.hpp index 60ac0747ab..564fc25eca 100644 --- a/modules/stitching/util.hpp +++ b/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()); } - int numVertices() const { return static_cast(edges_.size()); } - void addEdge(int from, int to, float weight); - - template - B forEach(B body) const; - - template - B walkBreadthFirst(int from, B body) const; + template B forEach(B body) const; + template B walkBreadthFirst(int from, B body) const; private: std::vector< std::list > 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 &corners, const std::vector &images); +cv::Rect resultRoi(const std::vector &corners, const std::vector &sizes); +cv::Point resultTl(const std::vector &corners); + + #include "util_inl.hpp" #endif // __OPENCV_STITCHING_UTIL_HPP__ diff --git a/modules/stitching/util_inl.hpp b/modules/stitching/util_inl.hpp index bdd1176bba..12f8396a88 100644 --- a/modules/stitching/util_inl.hpp +++ b/modules/stitching/util_inl.hpp @@ -105,11 +105,8 @@ double normL2sq(const cv::Mat &r) } -template -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__