From ace94d2ebf914544868fc267f7518042aaaab85c Mon Sep 17 00:00:00 2001 From: Alexey Spizhevoy Date: Thu, 9 Jun 2011 10:16:10 +0000 Subject: [PATCH] fixed bug in opencv_stitching (corrected resize images step), added matches checking (both 1->2 and 2->1 must be presented) --- modules/stitching/blenders.cpp | 618 ++++++++++----------- modules/stitching/main.cpp | 25 +- modules/stitching/matchers.cpp | 8 +- modules/stitching/motion_estimators.cpp | 950 ++++++++++++++++---------------- modules/stitching/seam_finders.cpp | 712 ++++++++++++------------ modules/stitching/warpers.hpp | 235 ++++---- modules/stitching/warpers_inl.hpp | 414 +++++++------- samples/gpu/driver_api_stereo_multi.cpp | 53 +- 8 files changed, 1532 insertions(+), 1483 deletions(-) diff --git a/modules/stitching/blenders.cpp b/modules/stitching/blenders.cpp index fe185a1..5797245 100644 --- a/modules/stitching/blenders.cpp +++ b/modules/stitching/blenders.cpp @@ -38,312 +38,312 @@ // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // -//M*/ -#include "blenders.hpp" -#include "util.hpp" - -using namespace std; -using namespace cv; - -static const float WEIGHT_EPS = 1e-5f; - -Ptr Blender::createDefault(int type) -{ - if (type == NO) - return new Blender(); - if (type == FEATHER) - return new FeatherBlender(); - if (type == MULTI_BAND) - return new MultiBandBlender(); - CV_Error(CV_StsBadArg, "unsupported blending method"); - return NULL; -} - - -void Blender::prepare(const vector &corners, const vector &sizes) -{ - prepare(resultRoi(corners, sizes)); -} - - -void Blender::prepare(Rect dst_roi) -{ - dst_.create(dst_roi.size(), CV_16SC3); - dst_.setTo(Scalar::all(0)); - dst_mask_.create(dst_roi.size(), CV_8U); - dst_mask_.setTo(Scalar::all(0)); - dst_roi_ = dst_roi; -} - - -void Blender::feed(const Mat &img, const Mat &mask, Point tl) -{ - CV_Assert(img.type() == CV_16SC3); - CV_Assert(mask.type() == CV_8U); - int dx = tl.x - dst_roi_.x; - int dy = tl.y - dst_roi_.y; - - for (int y = 0; y < img.rows; ++y) - { - const Point3_ *src_row = img.ptr >(y); - Point3_ *dst_row = dst_.ptr >(dy + y); - const uchar *mask_row = mask.ptr(y); - uchar *dst_mask_row = dst_mask_.ptr(dy + y); - - for (int x = 0; x < img.cols; ++x) - { - if (mask_row[x]) - dst_row[dx + x] = src_row[x]; - dst_mask_row[dx + x] |= mask_row[x]; - } - } -} - - -void Blender::blend(Mat &dst, Mat &dst_mask) -{ - dst_.setTo(Scalar::all(0), dst_mask_ == 0); - dst = dst_; - dst_mask = dst_mask_; - dst_.release(); - dst_mask_.release(); -} - - -void FeatherBlender::prepare(Rect dst_roi) -{ - Blender::prepare(dst_roi); - dst_weight_map_.create(dst_roi.size(), CV_32F); - dst_weight_map_.setTo(0); -} - - -void FeatherBlender::feed(const Mat &img, const Mat &mask, Point tl) -{ - CV_Assert(img.type() == CV_16SC3); - CV_Assert(mask.type() == CV_8U); - - createWeightMap(mask, sharpness_, weight_map_); - int dx = tl.x - dst_roi_.x; - int dy = tl.y - dst_roi_.y; - - for (int y = 0; y < img.rows; ++y) - { - const Point3_* src_row = img.ptr >(y); - Point3_* dst_row = dst_.ptr >(dy + y); - const float* weight_row = weight_map_.ptr(y); - float* dst_weight_row = dst_weight_map_.ptr(dy + y); - - for (int x = 0; x < img.cols; ++x) - { - dst_row[dx + x].x += static_cast(src_row[x].x * weight_row[x]); - dst_row[dx + x].y += static_cast(src_row[x].y * weight_row[x]); - dst_row[dx + x].z += static_cast(src_row[x].z * weight_row[x]); - dst_weight_row[dx + x] += weight_row[x]; - } - } -} - - -void FeatherBlender::blend(Mat &dst, Mat &dst_mask) -{ - normalize(dst_weight_map_, dst_); - dst_mask_ = dst_weight_map_ > WEIGHT_EPS; - Blender::blend(dst, dst_mask); -} - - -void MultiBandBlender::prepare(Rect dst_roi) -{ - dst_roi_final_ = dst_roi; - - // Crop unnecessary bands - double max_len = static_cast(max(dst_roi.width, dst_roi.height)); - num_bands_ = min(actual_num_bands_, static_cast(ceil(log(max_len) / log(2.0)))); - - // Add border to the final image, to ensure sizes are divided by (1 << num_bands_) - dst_roi.width += ((1 << num_bands_) - dst_roi.width % (1 << num_bands_)) % (1 << num_bands_); - dst_roi.height += ((1 << num_bands_) - dst_roi.height % (1 << num_bands_)) % (1 << num_bands_); - - Blender::prepare(dst_roi); - - dst_pyr_laplace_.resize(num_bands_ + 1); - dst_pyr_laplace_[0] = dst_; - - dst_band_weights_.resize(num_bands_ + 1); - dst_band_weights_[0].create(dst_roi.size(), CV_32F); - dst_band_weights_[0].setTo(0); - - for (int i = 1; i <= num_bands_; ++i) - { - dst_pyr_laplace_[i].create((dst_pyr_laplace_[i - 1].rows + 1) / 2, - (dst_pyr_laplace_[i - 1].cols + 1) / 2, CV_16SC3); - dst_band_weights_[i].create((dst_band_weights_[i - 1].rows + 1) / 2, - (dst_band_weights_[i - 1].cols + 1) / 2, CV_32F); - dst_pyr_laplace_[i].setTo(Scalar::all(0)); - dst_band_weights_[i].setTo(0); - } -} - - -void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl) -{ - CV_Assert(img.type() == CV_16SC3); - CV_Assert(mask.type() == CV_8U); - - // Keep source image in memory with small border - int gap = 3 * (1 << num_bands_); - Point tl_new(max(dst_roi_.x, tl.x - gap), - max(dst_roi_.y, tl.y - gap)); - Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap), - min(dst_roi_.br().y, tl.y + img.rows + gap)); - - // Ensure coordinates of top-left, bottom-right corners are divided by (1 << num_bands_). - // After that scale between layers is exactly 2. - // - // We do it to avoid interpolation problems when keeping sub-images only. There is no such problem when - // image is bordered to have size equal to the final image size, but this is too memory hungry approach. - tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_bands_); - tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_bands_); - int width = br_new.x - tl_new.x; - int height = br_new.y - tl_new.y; - width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_bands_); - height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_); - br_new.x = tl_new.x + width; - br_new.y = tl_new.y + height; - int dy = max(br_new.y - dst_roi_.br().y, 0); - int dx = max(br_new.x - dst_roi_.br().x, 0); - tl_new.x -= dx; br_new.x -= dx; - tl_new.y -= dy; br_new.y -= dy; - - int top = tl.y - tl_new.y; - int left = tl.x - tl_new.x; - int bottom = br_new.y - tl.y - img.rows; - int right = br_new.x - tl.x - img.cols; - - // Create the source image Laplacian pyramid - vector src_pyr_gauss(num_bands_ + 1); - copyMakeBorder(img, src_pyr_gauss[0], top, bottom, left, right, - BORDER_REFLECT); - for (int i = 0; i < num_bands_; ++i) - pyrDown(src_pyr_gauss[i], src_pyr_gauss[i + 1]); - vector src_pyr_laplace; - createLaplacePyr(src_pyr_gauss, src_pyr_laplace); - src_pyr_gauss.clear(); - - // Create the weight map Gaussian pyramid - Mat weight_map; - mask.convertTo(weight_map, CV_32F, 1./255.); - vector weight_pyr_gauss(num_bands_ + 1); - copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right, - BORDER_CONSTANT); - for (int i = 0; i < num_bands_; ++i) - pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]); - - int y_tl = tl_new.y - dst_roi_.y; - int y_br = br_new.y - dst_roi_.y; - int x_tl = tl_new.x - dst_roi_.x; - int x_br = br_new.x - dst_roi_.x; - - // Add weighted layer of the source image to the final Laplacian pyramid layer - for (int i = 0; i <= num_bands_; ++i) - { - for (int y = y_tl; y < y_br; ++y) - { - int y_ = y - y_tl; - const Point3_* src_row = src_pyr_laplace[i].ptr >(y_); - Point3_* dst_row = dst_pyr_laplace_[i].ptr >(y); - const float* weight_row = weight_pyr_gauss[i].ptr(y_); - float* dst_weight_row = dst_band_weights_[i].ptr(y); - - for (int x = x_tl; x < x_br; ++x) - { - int x_ = x - x_tl; - dst_row[x].x += static_cast(src_row[x_].x * weight_row[x_]); - dst_row[x].y += static_cast(src_row[x_].y * weight_row[x_]); - dst_row[x].z += static_cast(src_row[x_].z * weight_row[x_]); - dst_weight_row[x] += weight_row[x_]; - } - } - x_tl /= 2; y_tl /= 2; - x_br /= 2; y_br /= 2; - } -} - - -void MultiBandBlender::blend(Mat &dst, Mat &dst_mask) -{ - for (int i = 0; i <= num_bands_; ++i) - normalize(dst_band_weights_[i], dst_pyr_laplace_[i]); - - restoreImageFromLaplacePyr(dst_pyr_laplace_); - - dst_ = dst_pyr_laplace_[0]; - dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width)); - dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS; - dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width)); - dst_pyr_laplace_.clear(); - dst_band_weights_.clear(); - - Blender::blend(dst, dst_mask); -} - - -////////////////////////////////////////////////////////////////////////////// -// Auxiliary functions - -void normalize(const Mat& weight, Mat& src) -{ - CV_Assert(weight.type() == CV_32F); - CV_Assert(src.type() == CV_16SC3); - for (int y = 0; y < src.rows; ++y) - { - Point3_ *row = src.ptr >(y); - const float *weight_row = weight.ptr(y); - - for (int x = 0; x < src.cols; ++x) - { - row[x].x = static_cast(row[x].x / (weight_row[x] + WEIGHT_EPS)); - row[x].y = static_cast(row[x].y / (weight_row[x] + WEIGHT_EPS)); - row[x].z = static_cast(row[x].z / (weight_row[x] + WEIGHT_EPS)); - } - } -} - - -void createWeightMap(const Mat &mask, float sharpness, Mat &weight) -{ - CV_Assert(mask.type() == CV_8U); - distanceTransform(mask, weight, CV_DIST_L1, 3); - threshold(weight * sharpness, weight, 1.f, 1.f, THRESH_TRUNC); -} - - -void createLaplacePyr(const vector &pyr_gauss, vector &pyr_laplace) -{ - if (pyr_gauss.size() == 0) - return; - pyr_laplace.resize(pyr_gauss.size()); - Mat tmp; - for (size_t i = 0; i < pyr_laplace.size() - 1; ++i) - { - pyrUp(pyr_gauss[i + 1], tmp, pyr_gauss[i].size()); - subtract(pyr_gauss[i], tmp, pyr_laplace[i]); - } - pyr_laplace[pyr_laplace.size() - 1] = pyr_gauss[pyr_laplace.size() - 1].clone(); -} - - -void restoreImageFromLaplacePyr(vector &pyr) -{ - if (pyr.size() == 0) - return; - Mat tmp; - for (size_t i = pyr.size() - 1; i > 0; --i) - { - pyrUp(pyr[i], tmp, pyr[i - 1].size()); - add(tmp, pyr[i - 1], pyr[i - 1]); - } -} - - +//M*/ +#include "blenders.hpp" +#include "util.hpp" + +using namespace std; +using namespace cv; + +static const float WEIGHT_EPS = 1e-5f; + +Ptr Blender::createDefault(int type) +{ + if (type == NO) + return new Blender(); + if (type == FEATHER) + return new FeatherBlender(); + if (type == MULTI_BAND) + return new MultiBandBlender(); + CV_Error(CV_StsBadArg, "unsupported blending method"); + return NULL; +} + + +void Blender::prepare(const vector &corners, const vector &sizes) +{ + prepare(resultRoi(corners, sizes)); +} + + +void Blender::prepare(Rect dst_roi) +{ + dst_.create(dst_roi.size(), CV_16SC3); + dst_.setTo(Scalar::all(0)); + dst_mask_.create(dst_roi.size(), CV_8U); + dst_mask_.setTo(Scalar::all(0)); + dst_roi_ = dst_roi; +} + + +void Blender::feed(const Mat &img, const Mat &mask, Point tl) +{ + CV_Assert(img.type() == CV_16SC3); + CV_Assert(mask.type() == CV_8U); + int dx = tl.x - dst_roi_.x; + int dy = tl.y - dst_roi_.y; + + for (int y = 0; y < img.rows; ++y) + { + const Point3_ *src_row = img.ptr >(y); + Point3_ *dst_row = dst_.ptr >(dy + y); + const uchar *mask_row = mask.ptr(y); + uchar *dst_mask_row = dst_mask_.ptr(dy + y); + + for (int x = 0; x < img.cols; ++x) + { + if (mask_row[x]) + dst_row[dx + x] = src_row[x]; + dst_mask_row[dx + x] |= mask_row[x]; + } + } +} + + +void Blender::blend(Mat &dst, Mat &dst_mask) +{ + dst_.setTo(Scalar::all(0), dst_mask_ == 0); + dst = dst_; + dst_mask = dst_mask_; + dst_.release(); + dst_mask_.release(); +} + + +void FeatherBlender::prepare(Rect dst_roi) +{ + Blender::prepare(dst_roi); + dst_weight_map_.create(dst_roi.size(), CV_32F); + dst_weight_map_.setTo(0); +} + + +void FeatherBlender::feed(const Mat &img, const Mat &mask, Point tl) +{ + CV_Assert(img.type() == CV_16SC3); + CV_Assert(mask.type() == CV_8U); + + createWeightMap(mask, sharpness_, weight_map_); + int dx = tl.x - dst_roi_.x; + int dy = tl.y - dst_roi_.y; + + for (int y = 0; y < img.rows; ++y) + { + const Point3_* src_row = img.ptr >(y); + Point3_* dst_row = dst_.ptr >(dy + y); + const float* weight_row = weight_map_.ptr(y); + float* dst_weight_row = dst_weight_map_.ptr(dy + y); + + for (int x = 0; x < img.cols; ++x) + { + dst_row[dx + x].x += static_cast(src_row[x].x * weight_row[x]); + dst_row[dx + x].y += static_cast(src_row[x].y * weight_row[x]); + dst_row[dx + x].z += static_cast(src_row[x].z * weight_row[x]); + dst_weight_row[dx + x] += weight_row[x]; + } + } +} + + +void FeatherBlender::blend(Mat &dst, Mat &dst_mask) +{ + normalize(dst_weight_map_, dst_); + dst_mask_ = dst_weight_map_ > WEIGHT_EPS; + Blender::blend(dst, dst_mask); +} + + +void MultiBandBlender::prepare(Rect dst_roi) +{ + dst_roi_final_ = dst_roi; + + // Crop unnecessary bands + double max_len = static_cast(max(dst_roi.width, dst_roi.height)); + num_bands_ = min(actual_num_bands_, static_cast(ceil(log(max_len) / log(2.0)))); + + // Add border to the final image, to ensure sizes are divided by (1 << num_bands_) + dst_roi.width += ((1 << num_bands_) - dst_roi.width % (1 << num_bands_)) % (1 << num_bands_); + dst_roi.height += ((1 << num_bands_) - dst_roi.height % (1 << num_bands_)) % (1 << num_bands_); + + Blender::prepare(dst_roi); + + dst_pyr_laplace_.resize(num_bands_ + 1); + dst_pyr_laplace_[0] = dst_; + + dst_band_weights_.resize(num_bands_ + 1); + dst_band_weights_[0].create(dst_roi.size(), CV_32F); + dst_band_weights_[0].setTo(0); + + for (int i = 1; i <= num_bands_; ++i) + { + dst_pyr_laplace_[i].create((dst_pyr_laplace_[i - 1].rows + 1) / 2, + (dst_pyr_laplace_[i - 1].cols + 1) / 2, CV_16SC3); + dst_band_weights_[i].create((dst_band_weights_[i - 1].rows + 1) / 2, + (dst_band_weights_[i - 1].cols + 1) / 2, CV_32F); + dst_pyr_laplace_[i].setTo(Scalar::all(0)); + dst_band_weights_[i].setTo(0); + } +} + + +void MultiBandBlender::feed(const Mat &img, const Mat &mask, Point tl) +{ + CV_Assert(img.type() == CV_16SC3); + CV_Assert(mask.type() == CV_8U); + + // Keep source image in memory with small border + int gap = 3 * (1 << num_bands_); + Point tl_new(max(dst_roi_.x, tl.x - gap), + max(dst_roi_.y, tl.y - gap)); + Point br_new(min(dst_roi_.br().x, tl.x + img.cols + gap), + min(dst_roi_.br().y, tl.y + img.rows + gap)); + + // Ensure coordinates of top-left, bottom-right corners are divided by (1 << num_bands_). + // After that scale between layers is exactly 2. + // + // We do it to avoid interpolation problems when keeping sub-images only. There is no such problem when + // image is bordered to have size equal to the final image size, but this is too memory hungry approach. + tl_new.x = dst_roi_.x + (((tl_new.x - dst_roi_.x) >> num_bands_) << num_bands_); + tl_new.y = dst_roi_.y + (((tl_new.y - dst_roi_.y) >> num_bands_) << num_bands_); + int width = br_new.x - tl_new.x; + int height = br_new.y - tl_new.y; + width += ((1 << num_bands_) - width % (1 << num_bands_)) % (1 << num_bands_); + height += ((1 << num_bands_) - height % (1 << num_bands_)) % (1 << num_bands_); + br_new.x = tl_new.x + width; + br_new.y = tl_new.y + height; + int dy = max(br_new.y - dst_roi_.br().y, 0); + int dx = max(br_new.x - dst_roi_.br().x, 0); + tl_new.x -= dx; br_new.x -= dx; + tl_new.y -= dy; br_new.y -= dy; + + int top = tl.y - tl_new.y; + int left = tl.x - tl_new.x; + int bottom = br_new.y - tl.y - img.rows; + int right = br_new.x - tl.x - img.cols; + + // Create the source image Laplacian pyramid + vector src_pyr_gauss(num_bands_ + 1); + copyMakeBorder(img, src_pyr_gauss[0], top, bottom, left, right, + BORDER_REFLECT); + for (int i = 0; i < num_bands_; ++i) + pyrDown(src_pyr_gauss[i], src_pyr_gauss[i + 1]); + vector src_pyr_laplace; + createLaplacePyr(src_pyr_gauss, src_pyr_laplace); + src_pyr_gauss.clear(); + + // Create the weight map Gaussian pyramid + Mat weight_map; + mask.convertTo(weight_map, CV_32F, 1./255.); + vector weight_pyr_gauss(num_bands_ + 1); + copyMakeBorder(weight_map, weight_pyr_gauss[0], top, bottom, left, right, + BORDER_CONSTANT); + for (int i = 0; i < num_bands_; ++i) + pyrDown(weight_pyr_gauss[i], weight_pyr_gauss[i + 1]); + + int y_tl = tl_new.y - dst_roi_.y; + int y_br = br_new.y - dst_roi_.y; + int x_tl = tl_new.x - dst_roi_.x; + int x_br = br_new.x - dst_roi_.x; + + // Add weighted layer of the source image to the final Laplacian pyramid layer + for (int i = 0; i <= num_bands_; ++i) + { + for (int y = y_tl; y < y_br; ++y) + { + int y_ = y - y_tl; + const Point3_* src_row = src_pyr_laplace[i].ptr >(y_); + Point3_* dst_row = dst_pyr_laplace_[i].ptr >(y); + const float* weight_row = weight_pyr_gauss[i].ptr(y_); + float* dst_weight_row = dst_band_weights_[i].ptr(y); + + for (int x = x_tl; x < x_br; ++x) + { + int x_ = x - x_tl; + dst_row[x].x += static_cast(src_row[x_].x * weight_row[x_]); + dst_row[x].y += static_cast(src_row[x_].y * weight_row[x_]); + dst_row[x].z += static_cast(src_row[x_].z * weight_row[x_]); + dst_weight_row[x] += weight_row[x_]; + } + } + x_tl /= 2; y_tl /= 2; + x_br /= 2; y_br /= 2; + } +} + + +void MultiBandBlender::blend(Mat &dst, Mat &dst_mask) +{ + for (int i = 0; i <= num_bands_; ++i) + normalize(dst_band_weights_[i], dst_pyr_laplace_[i]); + + restoreImageFromLaplacePyr(dst_pyr_laplace_); + + dst_ = dst_pyr_laplace_[0]; + dst_ = dst_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width)); + dst_mask_ = dst_band_weights_[0] > WEIGHT_EPS; + dst_mask_ = dst_mask_(Range(0, dst_roi_final_.height), Range(0, dst_roi_final_.width)); + dst_pyr_laplace_.clear(); + dst_band_weights_.clear(); + + Blender::blend(dst, dst_mask); +} + + +////////////////////////////////////////////////////////////////////////////// +// Auxiliary functions + +void normalize(const Mat& weight, Mat& src) +{ + CV_Assert(weight.type() == CV_32F); + CV_Assert(src.type() == CV_16SC3); + for (int y = 0; y < src.rows; ++y) + { + Point3_ *row = src.ptr >(y); + const float *weight_row = weight.ptr(y); + + for (int x = 0; x < src.cols; ++x) + { + row[x].x = static_cast(row[x].x / (weight_row[x] + WEIGHT_EPS)); + row[x].y = static_cast(row[x].y / (weight_row[x] + WEIGHT_EPS)); + row[x].z = static_cast(row[x].z / (weight_row[x] + WEIGHT_EPS)); + } + } +} + + +void createWeightMap(const Mat &mask, float sharpness, Mat &weight) +{ + CV_Assert(mask.type() == CV_8U); + distanceTransform(mask, weight, CV_DIST_L1, 3); + threshold(weight * sharpness, weight, 1.f, 1.f, THRESH_TRUNC); +} + + +void createLaplacePyr(const vector &pyr_gauss, vector &pyr_laplace) +{ + if (pyr_gauss.size() == 0) + return; + pyr_laplace.resize(pyr_gauss.size()); + Mat tmp; + for (size_t i = 0; i < pyr_laplace.size() - 1; ++i) + { + pyrUp(pyr_gauss[i + 1], tmp, pyr_gauss[i].size()); + subtract(pyr_gauss[i], tmp, pyr_laplace[i]); + } + pyr_laplace[pyr_laplace.size() - 1] = pyr_gauss[pyr_laplace.size() - 1].clone(); +} + + +void restoreImageFromLaplacePyr(vector &pyr) +{ + if (pyr.size() == 0) + return; + Mat tmp; + for (size_t i = pyr.size() - 1; i > 0; --i) + { + pyrUp(pyr[i], tmp, pyr[i - 1].size()); + add(tmp, pyr[i - 1], pyr[i - 1]); + } +} + + diff --git a/modules/stitching/main.cpp b/modules/stitching/main.cpp index fe952e6..7c49228 100644 --- a/modules/stitching/main.cpp +++ b/modules/stitching/main.cpp @@ -75,7 +75,7 @@ void printUsage() " --work_megapix \n" " Resolution for image registration step. The default is 0.6 Mpx.\n" " --match_conf \n" - " Confidence for feature matching step. The default is 0.7.\n" + " Confidence for feature matching step. The default is 0.65.\n" " --conf_thresh \n" " Threshold for two images are from the same panorama confidence.\n" " The default is 1.0.\n" @@ -320,11 +320,14 @@ int main(int argc, char* argv[]) Mat full_img, img; vector images(num_images); + vector full_img_sizes(num_images); double seam_work_aspect = 1; for (int i = 0; i < num_images; ++i) { full_img = imread(img_names[i]); + full_img_sizes[i] = full_img.size(); + if (full_img.empty()) { LOGLN("Can't open image " << img_names[i]); @@ -376,14 +379,17 @@ int main(int argc, char* argv[]) vector indices = leaveBiggestComponent(features, pairwise_matches, conf_thresh); vector img_subset; vector img_names_subset; + vector full_img_sizes_subset; for (size_t i = 0; i < indices.size(); ++i) { img_names_subset.push_back(img_names[indices[i]]); img_subset.push_back(images[indices[i]]); + full_img_sizes_subset.push_back(full_img_sizes[indices[i]]); } images = img_subset; img_names = img_names_subset; + full_img_sizes = full_img_sizes_subset; // Check if we still have enough images num_images = static_cast(img_names.size()); @@ -519,16 +525,21 @@ int main(int argc, char* argv[]) 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)); + Size sz = full_img_sizes[i]; + if (abs(compose_scale - 1) > 1e-1) + { + sz.width = cvRound(full_img_sizes[i].width * compose_scale); + sz.height = cvRound(full_img_sizes[i].height * compose_scale); + } + Rect roi = warper->warpRoi(sz, static_cast(cameras[i].focal), cameras[i].R); + corners[i] = roi.tl(); + sizes[i] = roi.size(); } } if (abs(compose_scale - 1) > 1e-1) @@ -539,7 +550,7 @@ int main(int argc, char* argv[]) Size img_size = img.size(); // Warp the current image - warper->warp(img, static_cast(cameras[img_idx].focal), cameras[img_idx].R, + warper->warp(img, static_cast(cameras[img_idx].focal), cameras[img_idx].R, img_warped); // Warp the current image mask @@ -587,7 +598,7 @@ int main(int argc, char* argv[]) } Mat result, result_mask; - blender->blend(result, result_mask); + blender->blend(result, result_mask); LOGLN("Compositing, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec"); diff --git a/modules/stitching/matchers.cpp b/modules/stitching/matchers.cpp index 07c8fb9..69262eb 100644 --- a/modules/stitching/matchers.cpp +++ b/modules/stitching/matchers.cpp @@ -311,7 +311,7 @@ namespace const DMatch& m1 = pair_matches[i][1]; if (m0.distance < (1.f - match_conf_) * m1.distance) { - matches_info.matches.push_back(m0); + //matches_info.matches.push_back(m0); matches.insert(make_pair(m0.queryIdx, m0.trainIdx)); } } @@ -326,7 +326,7 @@ namespace const DMatch& m0 = pair_matches[i][0]; const DMatch& m1 = pair_matches[i][1]; if (m0.distance < (1.f - match_conf_) * m1.distance) - if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) == matches.end()) + if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) != matches.end()) matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance)); } } @@ -352,7 +352,7 @@ namespace const DMatch& m1 = pair_matches[i][1]; if (m0.distance < (1.f - match_conf_) * m1.distance) { - matches_info.matches.push_back(m0); + //matches_info.matches.push_back(m0); matches.insert(make_pair(m0.queryIdx, m0.trainIdx)); } } @@ -368,7 +368,7 @@ namespace const DMatch& m0 = pair_matches[i][0]; const DMatch& m1 = pair_matches[i][1]; if (m0.distance < (1.f - match_conf_) * m1.distance) - if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) == matches.end()) + if (matches.find(make_pair(m0.trainIdx, m0.queryIdx)) != matches.end()) matches_info.matches.push_back(DMatch(m0.trainIdx, m0.queryIdx, m0.distance)); } } diff --git a/modules/stitching/motion_estimators.cpp b/modules/stitching/motion_estimators.cpp index e1fac70..763c4a6 100644 --- a/modules/stitching/motion_estimators.cpp +++ b/modules/stitching/motion_estimators.cpp @@ -38,478 +38,478 @@ // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // -//M*/ -#include -#include "autocalib.hpp" -#include "motion_estimators.hpp" -#include "util.hpp" - -using namespace std; -using namespace cv; - - -////////////////////////////////////////////////////////////////////////////// - -CameraParams::CameraParams() : focal(1), R(Mat::eye(3, 3, CV_64F)), t(Mat::zeros(3, 1, CV_64F)) {} - -CameraParams::CameraParams(const CameraParams &other) { *this = other; } - -const CameraParams& CameraParams::operator =(const CameraParams &other) -{ - focal = other.focal; - R = other.R.clone(); - t = other.t.clone(); - return *this; -} - - -////////////////////////////////////////////////////////////////////////////// - -struct IncDistance -{ - IncDistance(vector &dists) : dists(&dists[0]) {} - void operator ()(const GraphEdge &edge) { dists[edge.to] = dists[edge.from] + 1; } - int* dists; -}; - - -struct CalcRotation -{ - CalcRotation(int num_images, const vector &pairwise_matches, vector &cameras) - : num_images(num_images), pairwise_matches(&pairwise_matches[0]), cameras(&cameras[0]) {} - - void operator ()(const GraphEdge &edge) - { - int pair_idx = edge.from * num_images + edge.to; - - double f_from = cameras[edge.from].focal; - double f_to = cameras[edge.to].focal; - - Mat K_from = Mat::eye(3, 3, CV_64F); - K_from.at(0, 0) = f_from; - K_from.at(1, 1) = f_from; - - Mat K_to = Mat::eye(3, 3, CV_64F); - K_to.at(0, 0) = f_to; - K_to.at(1, 1) = f_to; - - Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to; - cameras[edge.to].R = cameras[edge.from].R * R; - } - - int num_images; - const MatchesInfo* pairwise_matches; - CameraParams* cameras; -}; - - -void HomographyBasedEstimator::estimate(const vector &features, const vector &pairwise_matches, - vector &cameras) -{ - const int num_images = static_cast(features.size()); - - // Estimate focal length and set it for all cameras - vector focals; - estimateFocal(features, pairwise_matches, focals); - cameras.resize(num_images); - for (int i = 0; i < num_images; ++i) - cameras[i].focal = focals[i]; - - // Restore global motion - Graph span_tree; - vector span_tree_centers; - findMaxSpanningTree(num_images, pairwise_matches, span_tree, span_tree_centers); - span_tree.walkBreadthFirst(span_tree_centers[0], CalcRotation(num_images, pairwise_matches, cameras)); -} - - -////////////////////////////////////////////////////////////////////////////// - -void BundleAdjuster::estimate(const vector &features, const vector &pairwise_matches, - vector &cameras) -{ - num_images_ = static_cast(features.size()); - features_ = &features[0]; - pairwise_matches_ = &pairwise_matches[0]; - - // Prepare focals and rotations - cameras_.create(num_images_ * 4, 1, CV_64F); - SVD svd; - for (int i = 0; i < num_images_; ++i) - { - cameras_.at(i * 4, 0) = cameras[i].focal; - - svd(cameras[i].R, SVD::FULL_UV); - Mat R = svd.u * svd.vt; - if (determinant(R) < 0) - R *= -1; - - Mat rvec; - Rodrigues(R, rvec); CV_Assert(rvec.type() == CV_32F); - cameras_.at(i * 4 + 1, 0) = rvec.at(0, 0); - cameras_.at(i * 4 + 2, 0) = rvec.at(1, 0); - cameras_.at(i * 4 + 3, 0) = rvec.at(2, 0); - } - - // Select only consistent image pairs for futher adjustment - edges_.clear(); - for (int i = 0; i < num_images_ - 1; ++i) - { - for (int j = i + 1; j < num_images_; ++j) - { - const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; - if (matches_info.confidence > conf_thresh_) - edges_.push_back(make_pair(i, j)); - } - } - - // Compute number of correspondences - total_num_matches_ = 0; - for (size_t i = 0; i < edges_.size(); ++i) - total_num_matches_ += static_cast(pairwise_matches[edges_[i].first * num_images_ + edges_[i].second].num_inliers); - - CvLevMarq solver(num_images_ * 4, total_num_matches_ * 3, - cvTermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 1000, DBL_EPSILON)); - - CvMat matParams = cameras_; - cvCopy(&matParams, solver.param); - - int count = 0; - for(;;) - { - const CvMat* _param = 0; - CvMat* _J = 0; - CvMat* _err = 0; - - bool proceed = solver.update(_param, _J, _err); - - cvCopy( _param, &matParams ); - - if( !proceed || !_err ) - break; - - if( _J ) - { - calcJacobian(); - CvMat matJ = J_; - cvCopy( &matJ, _J ); - } - - if (_err) - { - calcError(err_); - LOG("."); - count++; - CvMat matErr = err_; - cvCopy( &matErr, _err ); - } - } - LOGLN(""); - LOGLN("Bundle adjustment, final error: " << sqrt(err_.dot(err_))); - LOGLN("Bundle adjustment, iterations done: " << count); - - // Obtain global motion - for (int i = 0; i < num_images_; ++i) - { - cameras[i].focal = cameras_.at(i * 4, 0); - Mat rvec(3, 1, CV_64F); - rvec.at(0, 0) = cameras_.at(i * 4 + 1, 0); - rvec.at(1, 0) = cameras_.at(i * 4 + 2, 0); - rvec.at(2, 0) = cameras_.at(i * 4 + 3, 0); - Rodrigues(rvec, cameras[i].R); - Mat Mf; - cameras[i].R.convertTo(Mf, CV_32F); - cameras[i].R = Mf; - } - - // Normalize motion to center image - Graph span_tree; - vector span_tree_centers; - findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_centers); - Mat R_inv = cameras[span_tree_centers[0]].R.inv(); - for (int i = 0; i < num_images_; ++i) - cameras[i].R = R_inv * cameras[i].R; -} - - -void BundleAdjuster::calcError(Mat &err) -{ - err.create(total_num_matches_ * 3, 1, CV_64F); - - int match_idx = 0; - for (size_t edge_idx = 0; edge_idx < edges_.size(); ++edge_idx) - { - int i = edges_[edge_idx].first; - int j = edges_[edge_idx].second; - double f1 = cameras_.at(i * 4, 0); - double f2 = cameras_.at(j * 4, 0); - double R1[9], R2[9]; - Mat R1_(3, 3, CV_64F, R1), R2_(3, 3, CV_64F, R2); - Mat rvec(3, 1, CV_64F); - rvec.at(0, 0) = cameras_.at(i * 4 + 1, 0); - rvec.at(1, 0) = cameras_.at(i * 4 + 2, 0); - rvec.at(2, 0) = cameras_.at(i * 4 + 3, 0); - Rodrigues(rvec, R1_); CV_Assert(R1_.type() == CV_64F); - rvec.at(0, 0) = cameras_.at(j * 4 + 1, 0); - rvec.at(1, 0) = cameras_.at(j * 4 + 2, 0); - rvec.at(2, 0) = cameras_.at(j * 4 + 3, 0); - Rodrigues(rvec, R2_); CV_Assert(R2_.type() == CV_64F); - - const ImageFeatures& features1 = features_[i]; - const ImageFeatures& features2 = features_[j]; - const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; - - for (size_t k = 0; k < matches_info.matches.size(); ++k) - { - if (!matches_info.inliers_mask[k]) - continue; - - const DMatch& m = matches_info.matches[k]; - - Point2d kp1 = features1.keypoints[m.queryIdx].pt; - kp1.x -= 0.5 * features1.img_size.width; - kp1.y -= 0.5 * features1.img_size.height; - Point2d kp2 = features2.keypoints[m.trainIdx].pt; - kp2.x -= 0.5 * features2.img_size.width; - kp2.y -= 0.5 * features2.img_size.height; - double len1 = sqrt(kp1.x * kp1.x + kp1.y * kp1.y + f1 * f1); - double len2 = sqrt(kp2.x * kp2.x + kp2.y * kp2.y + f2 * f2); - Point3d p1(kp1.x / len1, kp1.y / len1, f1 / len1); - Point3d p2(kp2.x / len2, kp2.y / len2, f2 / len2); - - Point3d d1(p1.x * R1[0] + p1.y * R1[1] + p1.z * R1[2], - p1.x * R1[3] + p1.y * R1[4] + p1.z * R1[5], - p1.x * R1[6] + p1.y * R1[7] + p1.z * R1[8]); - Point3d d2(p2.x * R2[0] + p2.y * R2[1] + p2.z * R2[2], - p2.x * R2[3] + p2.y * R2[4] + p2.z * R2[5], - p2.x * R2[6] + p2.y * R2[7] + p2.z * R2[8]); - - double mult = 1; - if (cost_space_ == FOCAL_RAY_SPACE) - mult = sqrt(f1 * f2); - err.at(3 * match_idx, 0) = mult * (d1.x - d2.x); - err.at(3 * match_idx + 1, 0) = mult * (d1.y - d2.y); - err.at(3 * match_idx + 2, 0) = mult * (d1.z - d2.z); - match_idx++; - } - } -} - - -void calcDeriv(const Mat &err1, const Mat &err2, double h, Mat res) -{ - for (int i = 0; i < err1.rows; ++i) - res.at(i, 0) = (err2.at(i, 0) - err1.at(i, 0)) / h; -} - - -void BundleAdjuster::calcJacobian() -{ - J_.create(total_num_matches_ * 3, num_images_ * 4, CV_64F); - - double f, r; - const double df = 0.001; // Focal length step - const double dr = 0.001; // Angle step - - for (int i = 0; i < num_images_; ++i) - { - f = cameras_.at(i * 4, 0); - cameras_.at(i * 4, 0) = f - df; - calcError(err1_); - cameras_.at(i * 4, 0) = f + df; - calcError(err2_); - calcDeriv(err1_, err2_, 2 * df, J_.col(i * 4)); - cameras_.at(i * 4, 0) = f; - - r = cameras_.at(i * 4 + 1, 0); - cameras_.at(i * 4 + 1, 0) = r - dr; - calcError(err1_); - cameras_.at(i * 4 + 1, 0) = r + dr; - calcError(err2_); - calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 1)); - cameras_.at(i * 4 + 1, 0) = r; - - r = cameras_.at(i * 4 + 2, 0); - cameras_.at(i * 4 + 2, 0) = r - dr; - calcError(err1_); - cameras_.at(i * 4 + 2, 0) = r + dr; - calcError(err2_); - calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 2)); - cameras_.at(i * 4 + 2, 0) = r; - - r = cameras_.at(i * 4 + 3, 0); - cameras_.at(i * 4 + 3, 0) = r - dr; - calcError(err1_); - cameras_.at(i * 4 + 3, 0) = r + dr; - calcError(err2_); - calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 3)); - cameras_.at(i * 4 + 3, 0) = r; - } -} - - -////////////////////////////////////////////////////////////////////////////// - -void waveCorrect(vector &rmats) -{ - float data[9]; - Mat r0(1, 3, CV_32F, data); - Mat r1(1, 3, CV_32F, data + 3); - Mat r2(1, 3, CV_32F, data + 6); - Mat R(3, 3, CV_32F, data); - - Mat cov = Mat::zeros(3, 3, CV_32F); - for (size_t i = 0; i < rmats.size(); ++i) - { - Mat r0 = rmats[i].col(0); - cov += r0 * r0.t(); - } - - SVD svd; - svd(cov, SVD::FULL_UV); - svd.vt.row(2).copyTo(r1); - if (determinant(svd.vt) < 0) r1 *= -1; - - Mat avgz = Mat::zeros(3, 1, CV_32F); - for (size_t i = 0; i < rmats.size(); ++i) - avgz += rmats[i].col(2); - r1.cross(avgz.t()).copyTo(r0); - normalize(r0, r0); - - r1.cross(r0).copyTo(r2); - if (determinant(R) < 0) R *= -1; - - for (size_t i = 0; i < rmats.size(); ++i) - rmats[i] = R * rmats[i]; -} - - -////////////////////////////////////////////////////////////////////////////// - -vector leaveBiggestComponent(vector &features, vector &pairwise_matches, - float conf_threshold) -{ - const int num_images = static_cast(features.size()); - - DjSets comps(num_images); - for (int i = 0; i < num_images; ++i) - { - for (int j = 0; j < num_images; ++j) - { - if (pairwise_matches[i*num_images + j].confidence < conf_threshold) - continue; - int comp1 = comps.find(i); - int comp2 = comps.find(j); - if (comp1 != comp2) - comps.merge(comp1, comp2); - } - } - - int max_comp = max_element(comps.size.begin(), comps.size.end()) - comps.size.begin(); - - vector indices; - vector indices_removed; - for (int i = 0; i < num_images; ++i) - if (comps.find(i) == max_comp) - indices.push_back(i); - else - indices_removed.push_back(i); - - vector features_subset; - vector pairwise_matches_subset; - for (size_t i = 0; i < indices.size(); ++i) - { - features_subset.push_back(features[indices[i]]); - for (size_t j = 0; j < indices.size(); ++j) - { - pairwise_matches_subset.push_back(pairwise_matches[indices[i]*num_images + indices[j]]); - pairwise_matches_subset.back().src_img_idx = i; - pairwise_matches_subset.back().dst_img_idx = j; - } - } - - if (static_cast(features_subset.size()) == num_images) - return indices; - - LOG("Removed some images, because can't match them: ("); - LOG(indices_removed[0]+1); - for (size_t i = 1; i < indices_removed.size(); ++i) - LOG(", " << indices_removed[i]+1); - LOGLN(")"); - - features = features_subset; - pairwise_matches = pairwise_matches_subset; - - return indices; -} - - -void findMaxSpanningTree(int num_images, const vector &pairwise_matches, - Graph &span_tree, vector ¢ers) -{ - Graph graph(num_images); - vector edges; - - // Construct images graph and remember its edges - for (int i = 0; i < num_images; ++i) - { - for (int j = 0; j < num_images; ++j) - { - if (pairwise_matches[i * num_images + j].H.empty()) - continue; - float conf = static_cast(pairwise_matches[i * num_images + j].num_inliers); - graph.addEdge(i, j, conf); - edges.push_back(GraphEdge(i, j, conf)); - } - } - - DjSets comps(num_images); - span_tree.create(num_images); - vector span_tree_powers(num_images, 0); - - // Find maximum spanning tree - sort(edges.begin(), edges.end(), greater()); - for (size_t i = 0; i < edges.size(); ++i) - { - int comp1 = comps.find(edges[i].from); - int comp2 = comps.find(edges[i].to); - if (comp1 != comp2) - { - comps.merge(comp1, comp2); - span_tree.addEdge(edges[i].from, edges[i].to, edges[i].weight); - span_tree.addEdge(edges[i].to, edges[i].from, edges[i].weight); - span_tree_powers[edges[i].from]++; - span_tree_powers[edges[i].to]++; - } - } - - // Find spanning tree leafs - vector span_tree_leafs; - for (int i = 0; i < num_images; ++i) - if (span_tree_powers[i] == 1) - span_tree_leafs.push_back(i); - - // Find maximum distance from each spanning tree vertex - vector max_dists(num_images, 0); - vector cur_dists; - for (size_t i = 0; i < span_tree_leafs.size(); ++i) - { - cur_dists.assign(num_images, 0); - span_tree.walkBreadthFirst(span_tree_leafs[i], IncDistance(cur_dists)); - for (int j = 0; j < num_images; ++j) - max_dists[j] = max(max_dists[j], cur_dists[j]); - } - - // Find min-max distance - int min_max_dist = max_dists[0]; - for (int i = 1; i < num_images; ++i) - if (min_max_dist > max_dists[i]) - min_max_dist = max_dists[i]; - - // Find spanning tree centers - centers.clear(); - for (int i = 0; i < num_images; ++i) - if (max_dists[i] == min_max_dist) - centers.push_back(i); - CV_Assert(centers.size() > 0 && centers.size() <= 2); -} +//M*/ +#include +#include "autocalib.hpp" +#include "motion_estimators.hpp" +#include "util.hpp" + +using namespace std; +using namespace cv; + + +////////////////////////////////////////////////////////////////////////////// + +CameraParams::CameraParams() : focal(1), R(Mat::eye(3, 3, CV_64F)), t(Mat::zeros(3, 1, CV_64F)) {} + +CameraParams::CameraParams(const CameraParams &other) { *this = other; } + +const CameraParams& CameraParams::operator =(const CameraParams &other) +{ + focal = other.focal; + R = other.R.clone(); + t = other.t.clone(); + return *this; +} + + +////////////////////////////////////////////////////////////////////////////// + +struct IncDistance +{ + IncDistance(vector &dists) : dists(&dists[0]) {} + void operator ()(const GraphEdge &edge) { dists[edge.to] = dists[edge.from] + 1; } + int* dists; +}; + + +struct CalcRotation +{ + CalcRotation(int num_images, const vector &pairwise_matches, vector &cameras) + : num_images(num_images), pairwise_matches(&pairwise_matches[0]), cameras(&cameras[0]) {} + + void operator ()(const GraphEdge &edge) + { + int pair_idx = edge.from * num_images + edge.to; + + double f_from = cameras[edge.from].focal; + double f_to = cameras[edge.to].focal; + + Mat K_from = Mat::eye(3, 3, CV_64F); + K_from.at(0, 0) = f_from; + K_from.at(1, 1) = f_from; + + Mat K_to = Mat::eye(3, 3, CV_64F); + K_to.at(0, 0) = f_to; + K_to.at(1, 1) = f_to; + + Mat R = K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to; + cameras[edge.to].R = cameras[edge.from].R * R; + } + + int num_images; + const MatchesInfo* pairwise_matches; + CameraParams* cameras; +}; + + +void HomographyBasedEstimator::estimate(const vector &features, const vector &pairwise_matches, + vector &cameras) +{ + const int num_images = static_cast(features.size()); + + // Estimate focal length and set it for all cameras + vector focals; + estimateFocal(features, pairwise_matches, focals); + cameras.resize(num_images); + for (int i = 0; i < num_images; ++i) + cameras[i].focal = focals[i]; + + // Restore global motion + Graph span_tree; + vector span_tree_centers; + findMaxSpanningTree(num_images, pairwise_matches, span_tree, span_tree_centers); + span_tree.walkBreadthFirst(span_tree_centers[0], CalcRotation(num_images, pairwise_matches, cameras)); +} + + +////////////////////////////////////////////////////////////////////////////// + +void BundleAdjuster::estimate(const vector &features, const vector &pairwise_matches, + vector &cameras) +{ + num_images_ = static_cast(features.size()); + features_ = &features[0]; + pairwise_matches_ = &pairwise_matches[0]; + + // Prepare focals and rotations + cameras_.create(num_images_ * 4, 1, CV_64F); + SVD svd; + for (int i = 0; i < num_images_; ++i) + { + cameras_.at(i * 4, 0) = cameras[i].focal; + + svd(cameras[i].R, SVD::FULL_UV); + Mat R = svd.u * svd.vt; + if (determinant(R) < 0) + R *= -1; + + Mat rvec; + Rodrigues(R, rvec); CV_Assert(rvec.type() == CV_32F); + cameras_.at(i * 4 + 1, 0) = rvec.at(0, 0); + cameras_.at(i * 4 + 2, 0) = rvec.at(1, 0); + cameras_.at(i * 4 + 3, 0) = rvec.at(2, 0); + } + + // Select only consistent image pairs for futher adjustment + edges_.clear(); + for (int i = 0; i < num_images_ - 1; ++i) + { + for (int j = i + 1; j < num_images_; ++j) + { + const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; + if (matches_info.confidence > conf_thresh_) + edges_.push_back(make_pair(i, j)); + } + } + + // Compute number of correspondences + total_num_matches_ = 0; + for (size_t i = 0; i < edges_.size(); ++i) + total_num_matches_ += static_cast(pairwise_matches[edges_[i].first * num_images_ + edges_[i].second].num_inliers); + + CvLevMarq solver(num_images_ * 4, total_num_matches_ * 3, + cvTermCriteria(CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 1000, DBL_EPSILON)); + + CvMat matParams = cameras_; + cvCopy(&matParams, solver.param); + + int count = 0; + for(;;) + { + const CvMat* _param = 0; + CvMat* _J = 0; + CvMat* _err = 0; + + bool proceed = solver.update(_param, _J, _err); + + cvCopy( _param, &matParams ); + + if( !proceed || !_err ) + break; + + if( _J ) + { + calcJacobian(); + CvMat matJ = J_; + cvCopy( &matJ, _J ); + } + + if (_err) + { + calcError(err_); + LOG("."); + count++; + CvMat matErr = err_; + cvCopy( &matErr, _err ); + } + } + LOGLN(""); + LOGLN("Bundle adjustment, final error: " << sqrt(err_.dot(err_))); + LOGLN("Bundle adjustment, iterations done: " << count); + + // Obtain global motion + for (int i = 0; i < num_images_; ++i) + { + cameras[i].focal = cameras_.at(i * 4, 0); + Mat rvec(3, 1, CV_64F); + rvec.at(0, 0) = cameras_.at(i * 4 + 1, 0); + rvec.at(1, 0) = cameras_.at(i * 4 + 2, 0); + rvec.at(2, 0) = cameras_.at(i * 4 + 3, 0); + Rodrigues(rvec, cameras[i].R); + Mat Mf; + cameras[i].R.convertTo(Mf, CV_32F); + cameras[i].R = Mf; + } + + // Normalize motion to center image + Graph span_tree; + vector span_tree_centers; + findMaxSpanningTree(num_images_, pairwise_matches, span_tree, span_tree_centers); + Mat R_inv = cameras[span_tree_centers[0]].R.inv(); + for (int i = 0; i < num_images_; ++i) + cameras[i].R = R_inv * cameras[i].R; +} + + +void BundleAdjuster::calcError(Mat &err) +{ + err.create(total_num_matches_ * 3, 1, CV_64F); + + int match_idx = 0; + for (size_t edge_idx = 0; edge_idx < edges_.size(); ++edge_idx) + { + int i = edges_[edge_idx].first; + int j = edges_[edge_idx].second; + double f1 = cameras_.at(i * 4, 0); + double f2 = cameras_.at(j * 4, 0); + double R1[9], R2[9]; + Mat R1_(3, 3, CV_64F, R1), R2_(3, 3, CV_64F, R2); + Mat rvec(3, 1, CV_64F); + rvec.at(0, 0) = cameras_.at(i * 4 + 1, 0); + rvec.at(1, 0) = cameras_.at(i * 4 + 2, 0); + rvec.at(2, 0) = cameras_.at(i * 4 + 3, 0); + Rodrigues(rvec, R1_); CV_Assert(R1_.type() == CV_64F); + rvec.at(0, 0) = cameras_.at(j * 4 + 1, 0); + rvec.at(1, 0) = cameras_.at(j * 4 + 2, 0); + rvec.at(2, 0) = cameras_.at(j * 4 + 3, 0); + Rodrigues(rvec, R2_); CV_Assert(R2_.type() == CV_64F); + + const ImageFeatures& features1 = features_[i]; + const ImageFeatures& features2 = features_[j]; + const MatchesInfo& matches_info = pairwise_matches_[i * num_images_ + j]; + + for (size_t k = 0; k < matches_info.matches.size(); ++k) + { + if (!matches_info.inliers_mask[k]) + continue; + + const DMatch& m = matches_info.matches[k]; + + Point2d kp1 = features1.keypoints[m.queryIdx].pt; + kp1.x -= 0.5 * features1.img_size.width; + kp1.y -= 0.5 * features1.img_size.height; + Point2d kp2 = features2.keypoints[m.trainIdx].pt; + kp2.x -= 0.5 * features2.img_size.width; + kp2.y -= 0.5 * features2.img_size.height; + double len1 = sqrt(kp1.x * kp1.x + kp1.y * kp1.y + f1 * f1); + double len2 = sqrt(kp2.x * kp2.x + kp2.y * kp2.y + f2 * f2); + Point3d p1(kp1.x / len1, kp1.y / len1, f1 / len1); + Point3d p2(kp2.x / len2, kp2.y / len2, f2 / len2); + + Point3d d1(p1.x * R1[0] + p1.y * R1[1] + p1.z * R1[2], + p1.x * R1[3] + p1.y * R1[4] + p1.z * R1[5], + p1.x * R1[6] + p1.y * R1[7] + p1.z * R1[8]); + Point3d d2(p2.x * R2[0] + p2.y * R2[1] + p2.z * R2[2], + p2.x * R2[3] + p2.y * R2[4] + p2.z * R2[5], + p2.x * R2[6] + p2.y * R2[7] + p2.z * R2[8]); + + double mult = 1; + if (cost_space_ == FOCAL_RAY_SPACE) + mult = sqrt(f1 * f2); + err.at(3 * match_idx, 0) = mult * (d1.x - d2.x); + err.at(3 * match_idx + 1, 0) = mult * (d1.y - d2.y); + err.at(3 * match_idx + 2, 0) = mult * (d1.z - d2.z); + match_idx++; + } + } +} + + +void calcDeriv(const Mat &err1, const Mat &err2, double h, Mat res) +{ + for (int i = 0; i < err1.rows; ++i) + res.at(i, 0) = (err2.at(i, 0) - err1.at(i, 0)) / h; +} + + +void BundleAdjuster::calcJacobian() +{ + J_.create(total_num_matches_ * 3, num_images_ * 4, CV_64F); + + double f, r; + const double df = 0.001; // Focal length step + const double dr = 0.001; // Angle step + + for (int i = 0; i < num_images_; ++i) + { + f = cameras_.at(i * 4, 0); + cameras_.at(i * 4, 0) = f - df; + calcError(err1_); + cameras_.at(i * 4, 0) = f + df; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * df, J_.col(i * 4)); + cameras_.at(i * 4, 0) = f; + + r = cameras_.at(i * 4 + 1, 0); + cameras_.at(i * 4 + 1, 0) = r - dr; + calcError(err1_); + cameras_.at(i * 4 + 1, 0) = r + dr; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 1)); + cameras_.at(i * 4 + 1, 0) = r; + + r = cameras_.at(i * 4 + 2, 0); + cameras_.at(i * 4 + 2, 0) = r - dr; + calcError(err1_); + cameras_.at(i * 4 + 2, 0) = r + dr; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 2)); + cameras_.at(i * 4 + 2, 0) = r; + + r = cameras_.at(i * 4 + 3, 0); + cameras_.at(i * 4 + 3, 0) = r - dr; + calcError(err1_); + cameras_.at(i * 4 + 3, 0) = r + dr; + calcError(err2_); + calcDeriv(err1_, err2_, 2 * dr, J_.col(i * 4 + 3)); + cameras_.at(i * 4 + 3, 0) = r; + } +} + + +////////////////////////////////////////////////////////////////////////////// + +void waveCorrect(vector &rmats) +{ + float data[9]; + Mat r0(1, 3, CV_32F, data); + Mat r1(1, 3, CV_32F, data + 3); + Mat r2(1, 3, CV_32F, data + 6); + Mat R(3, 3, CV_32F, data); + + Mat cov = Mat::zeros(3, 3, CV_32F); + for (size_t i = 0; i < rmats.size(); ++i) + { + Mat r0 = rmats[i].col(0); + cov += r0 * r0.t(); + } + + SVD svd; + svd(cov, SVD::FULL_UV); + svd.vt.row(2).copyTo(r1); + if (determinant(svd.vt) < 0) r1 *= -1; + + Mat avgz = Mat::zeros(3, 1, CV_32F); + for (size_t i = 0; i < rmats.size(); ++i) + avgz += rmats[i].col(2); + r1.cross(avgz.t()).copyTo(r0); + normalize(r0, r0); + + r1.cross(r0).copyTo(r2); + if (determinant(R) < 0) R *= -1; + + for (size_t i = 0; i < rmats.size(); ++i) + rmats[i] = R * rmats[i]; +} + + +////////////////////////////////////////////////////////////////////////////// + +vector leaveBiggestComponent(vector &features, vector &pairwise_matches, + float conf_threshold) +{ + const int num_images = static_cast(features.size()); + + DjSets comps(num_images); + for (int i = 0; i < num_images; ++i) + { + for (int j = 0; j < num_images; ++j) + { + if (pairwise_matches[i*num_images + j].confidence < conf_threshold) + continue; + int comp1 = comps.find(i); + int comp2 = comps.find(j); + if (comp1 != comp2) + comps.merge(comp1, comp2); + } + } + + int max_comp = max_element(comps.size.begin(), comps.size.end()) - comps.size.begin(); + + vector indices; + vector indices_removed; + for (int i = 0; i < num_images; ++i) + if (comps.find(i) == max_comp) + indices.push_back(i); + else + indices_removed.push_back(i); + + vector features_subset; + vector pairwise_matches_subset; + for (size_t i = 0; i < indices.size(); ++i) + { + features_subset.push_back(features[indices[i]]); + for (size_t j = 0; j < indices.size(); ++j) + { + pairwise_matches_subset.push_back(pairwise_matches[indices[i]*num_images + indices[j]]); + pairwise_matches_subset.back().src_img_idx = i; + pairwise_matches_subset.back().dst_img_idx = j; + } + } + + if (static_cast(features_subset.size()) == num_images) + return indices; + + LOG("Removed some images, because can't match them: ("); + LOG(indices_removed[0]+1); + for (size_t i = 1; i < indices_removed.size(); ++i) + LOG(", " << indices_removed[i]+1); + LOGLN("). Try decrease --match_conf value."); + + features = features_subset; + pairwise_matches = pairwise_matches_subset; + + return indices; +} + + +void findMaxSpanningTree(int num_images, const vector &pairwise_matches, + Graph &span_tree, vector ¢ers) +{ + Graph graph(num_images); + vector edges; + + // Construct images graph and remember its edges + for (int i = 0; i < num_images; ++i) + { + for (int j = 0; j < num_images; ++j) + { + if (pairwise_matches[i * num_images + j].H.empty()) + continue; + float conf = static_cast(pairwise_matches[i * num_images + j].num_inliers); + graph.addEdge(i, j, conf); + edges.push_back(GraphEdge(i, j, conf)); + } + } + + DjSets comps(num_images); + span_tree.create(num_images); + vector span_tree_powers(num_images, 0); + + // Find maximum spanning tree + sort(edges.begin(), edges.end(), greater()); + for (size_t i = 0; i < edges.size(); ++i) + { + int comp1 = comps.find(edges[i].from); + int comp2 = comps.find(edges[i].to); + if (comp1 != comp2) + { + comps.merge(comp1, comp2); + span_tree.addEdge(edges[i].from, edges[i].to, edges[i].weight); + span_tree.addEdge(edges[i].to, edges[i].from, edges[i].weight); + span_tree_powers[edges[i].from]++; + span_tree_powers[edges[i].to]++; + } + } + + // Find spanning tree leafs + vector span_tree_leafs; + for (int i = 0; i < num_images; ++i) + if (span_tree_powers[i] == 1) + span_tree_leafs.push_back(i); + + // Find maximum distance from each spanning tree vertex + vector max_dists(num_images, 0); + vector cur_dists; + for (size_t i = 0; i < span_tree_leafs.size(); ++i) + { + cur_dists.assign(num_images, 0); + span_tree.walkBreadthFirst(span_tree_leafs[i], IncDistance(cur_dists)); + for (int j = 0; j < num_images; ++j) + max_dists[j] = max(max_dists[j], cur_dists[j]); + } + + // Find min-max distance + int min_max_dist = max_dists[0]; + for (int i = 1; i < num_images; ++i) + if (min_max_dist > max_dists[i]) + min_max_dist = max_dists[i]; + + // Find spanning tree centers + centers.clear(); + for (int i = 0; i < num_images; ++i) + if (max_dists[i] == min_max_dist) + centers.push_back(i); + CV_Assert(centers.size() > 0 && centers.size() <= 2); +} diff --git a/modules/stitching/seam_finders.cpp b/modules/stitching/seam_finders.cpp index c2eb8eb..efcd86e 100644 --- a/modules/stitching/seam_finders.cpp +++ b/modules/stitching/seam_finders.cpp @@ -38,365 +38,365 @@ // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // -//M*/ -#include "seam_finders.hpp" -#include "util.hpp" - -using namespace std; -using namespace cv; - - -Ptr SeamFinder::createDefault(int type) -{ - if (type == NO) - return new NoSeamFinder(); - if (type == VORONOI) - return new VoronoiSeamFinder(); - if (type == GC_COLOR) - return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR); - if (type == GC_COLOR_GRAD) - return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR_GRAD); - CV_Error(CV_StsBadArg, "unsupported seam finding method"); - return NULL; -} - - -void PairwiseSeamFinder::find(const vector &src, const vector &corners, - vector &masks) -{ - if (src.size() == 0) - return; - - images_ = src; - corners_ = corners; - masks_ = masks; - - for (size_t i = 0; i < src.size() - 1; ++i) - { - for (size_t j = i + 1; j < src.size(); ++j) - { - Rect roi; - if (overlapRoi(corners[i], corners[j], src[i].size(), src[j].size(), roi)) - findInPair(i, j, roi); - } - } -} - - -void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi) -{ - 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); - - Mat img1 = images_[first], img2 = images_[second]; - Mat mask1 = masks_[first], mask2 = masks_[second]; - Point tl1 = corners_[first], tl2 = corners_[second]; - - // Cut submasks with some gap - for (int y = -gap; y < roi.height + gap; ++y) - { - for (int x = -gap; x < roi.width + gap; ++x) - { - int y1 = roi.y - tl1.y + y; - int x1 = roi.x - tl1.x + x; - if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) - submask1.at(y + gap, x + gap) = mask1.at(y1, x1); - else - submask1.at(y + gap, x + gap) = 0; - - int y2 = roi.y - tl2.y + y; - int x2 = roi.x - tl2.x + x; - if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) - submask2.at(y + gap, x + gap) = mask2.at(y2, x2); - else - submask2.at(y + gap, x + gap) = 0; - } - } - - Mat collision = (submask1 != 0) & (submask2 != 0); - Mat unique1 = submask1.clone(); unique1.setTo(0, collision); - Mat unique2 = submask2.clone(); unique2.setTo(0, collision); - - Mat dist1, dist2; - distanceTransform(unique1 == 0, dist1, CV_DIST_L1, 3); - distanceTransform(unique2 == 0, dist2, CV_DIST_L1, 3); - - Mat seam = dist1 < dist2; - - for (int y = 0; y < roi.height; ++y) - { - for (int x = 0; x < roi.width; ++x) - { - if (seam.at(y + gap, x + gap)) - mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; - else - mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; - } - } -} - - -class GraphCutSeamFinder::Impl : public PairwiseSeamFinder -{ -public: - Impl(int cost_type, float terminal_cost, float bad_region_penalty) - : cost_type_(cost_type), terminal_cost_(terminal_cost), bad_region_penalty_(bad_region_penalty) {} - - void find(const vector &src, const vector &corners, vector &masks); - void findInPair(size_t first, size_t second, Rect roi); - -private: +//M*/ +#include "seam_finders.hpp" +#include "util.hpp" + +using namespace std; +using namespace cv; + + +Ptr SeamFinder::createDefault(int type) +{ + if (type == NO) + return new NoSeamFinder(); + if (type == VORONOI) + return new VoronoiSeamFinder(); + if (type == GC_COLOR) + return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR); + if (type == GC_COLOR_GRAD) + return new GraphCutSeamFinder(GraphCutSeamFinder::COST_COLOR_GRAD); + CV_Error(CV_StsBadArg, "unsupported seam finding method"); + return NULL; +} + + +void PairwiseSeamFinder::find(const vector &src, const vector &corners, + vector &masks) +{ + if (src.size() == 0) + return; + + images_ = src; + corners_ = corners; + masks_ = masks; + + for (size_t i = 0; i < src.size() - 1; ++i) + { + for (size_t j = i + 1; j < src.size(); ++j) + { + Rect roi; + if (overlapRoi(corners[i], corners[j], src[i].size(), src[j].size(), roi)) + findInPair(i, j, roi); + } + } +} + + +void VoronoiSeamFinder::findInPair(size_t first, size_t second, Rect roi) +{ + 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); + + Mat img1 = images_[first], img2 = images_[second]; + Mat mask1 = masks_[first], mask2 = masks_[second]; + Point tl1 = corners_[first], tl2 = corners_[second]; + + // Cut submasks with some gap + for (int y = -gap; y < roi.height + gap; ++y) + { + for (int x = -gap; x < roi.width + gap; ++x) + { + int y1 = roi.y - tl1.y + y; + int x1 = roi.x - tl1.x + x; + if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) + submask1.at(y + gap, x + gap) = mask1.at(y1, x1); + else + submask1.at(y + gap, x + gap) = 0; + + int y2 = roi.y - tl2.y + y; + int x2 = roi.x - tl2.x + x; + if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) + submask2.at(y + gap, x + gap) = mask2.at(y2, x2); + else + submask2.at(y + gap, x + gap) = 0; + } + } + + Mat collision = (submask1 != 0) & (submask2 != 0); + Mat unique1 = submask1.clone(); unique1.setTo(0, collision); + Mat unique2 = submask2.clone(); unique2.setTo(0, collision); + + Mat dist1, dist2; + distanceTransform(unique1 == 0, dist1, CV_DIST_L1, 3); + distanceTransform(unique2 == 0, dist2, CV_DIST_L1, 3); + + Mat seam = dist1 < dist2; + + for (int y = 0; y < roi.height; ++y) + { + for (int x = 0; x < roi.width; ++x) + { + if (seam.at(y + gap, x + gap)) + mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; + else + mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; + } + } +} + + +class GraphCutSeamFinder::Impl : public PairwiseSeamFinder +{ +public: + Impl(int cost_type, float terminal_cost, float bad_region_penalty) + : cost_type_(cost_type), terminal_cost_(terminal_cost), bad_region_penalty_(bad_region_penalty) {} + + void find(const vector &src, const vector &corners, vector &masks); + void findInPair(size_t first, size_t second, Rect roi); + +private: void setGraphWeightsColor(const Mat &img1, const Mat &img2, - const Mat &mask1, const Mat &mask2, GCGraph &graph); + const Mat &mask1, const Mat &mask2, GCGraph &graph); void setGraphWeightsColorGrad(const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, - const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, - GCGraph &graph); - - vector dx_, dy_; - int cost_type_; - float terminal_cost_; - float bad_region_penalty_; -}; - - -void GraphCutSeamFinder::Impl::find(const vector &src, const vector &corners, - vector &masks) -{ - // Compute gradients - dx_.resize(src.size()); - dy_.resize(src.size()); - Mat dx, dy; - for (size_t i = 0; i < src.size(); ++i) - { - CV_Assert(src[i].channels() == 3); - Sobel(src[i], dx, CV_32F, 1, 0); - Sobel(src[i], dy, CV_32F, 0, 1); - dx_[i].create(src[i].size(), CV_32F); - dy_[i].create(src[i].size(), CV_32F); - for (int y = 0; y < src[i].rows; ++y) - { - const Point3f* dx_row = dx.ptr(y); - const Point3f* dy_row = dy.ptr(y); - float* dx_row_ = dx_[i].ptr(y); - float* dy_row_ = dy_[i].ptr(y); - for (int x = 0; x < src[i].cols; ++x) - { - dx_row_[x] = normL2(dx_row[x]); - dy_row_[x] = normL2(dy_row[x]); - } - } - } - PairwiseSeamFinder::find(src, corners, masks); -} - - + const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, + GCGraph &graph); + + vector dx_, dy_; + int cost_type_; + float terminal_cost_; + float bad_region_penalty_; +}; + + +void GraphCutSeamFinder::Impl::find(const vector &src, const vector &corners, + vector &masks) +{ + // Compute gradients + dx_.resize(src.size()); + dy_.resize(src.size()); + Mat dx, dy; + for (size_t i = 0; i < src.size(); ++i) + { + CV_Assert(src[i].channels() == 3); + Sobel(src[i], dx, CV_32F, 1, 0); + Sobel(src[i], dy, CV_32F, 0, 1); + dx_[i].create(src[i].size(), CV_32F); + dy_[i].create(src[i].size(), CV_32F); + for (int y = 0; y < src[i].rows; ++y) + { + const Point3f* dx_row = dx.ptr(y); + const Point3f* dy_row = dy.ptr(y); + float* dx_row_ = dx_[i].ptr(y); + float* dy_row_ = dy_[i].ptr(y); + for (int x = 0; x < src[i].cols; ++x) + { + dx_row_[x] = normL2(dx_row[x]); + dy_row_[x] = normL2(dy_row[x]); + } + } + } + PairwiseSeamFinder::find(src, corners, masks); +} + + void GraphCutSeamFinder::Impl::setGraphWeightsColor(const Mat &img1, const Mat &img2, - const Mat &mask1, const Mat &mask2, GCGraph &graph) -{ - const Size img_size = img1.size(); - - // Set terminal weights - for (int y = 0; y < img_size.height; ++y) - { - for (int x = 0; x < img_size.width; ++x) - { - int v = graph.addVtx(); - graph.addTermWeights(v, mask1.at(y, x) ? terminal_cost_ : 0.f, - mask2.at(y, x) ? terminal_cost_ : 0.f); - } - } - - // Set regular edge weights - const float weight_eps = 1e-3f; - for (int y = 0; y < img_size.height; ++y) - { - for (int x = 0; x < img_size.width; ++x) - { - int v = y * img_size.width + x; - if (x < img_size.width - 1) - { - float weight = normL2(img1.at(y, x), img2.at(y, x)) + - normL2(img1.at(y, x + 1), img2.at(y, x + 1)) + - weight_eps; - if (!mask1.at(y, x) || !mask1.at(y, x + 1) || - !mask2.at(y, x) || !mask2.at(y, x + 1)) - weight += bad_region_penalty_; - graph.addEdges(v, v + 1, weight, weight); - } - if (y < img_size.height - 1) - { - float weight = normL2(img1.at(y, x), img2.at(y, x)) + - normL2(img1.at(y + 1, x), img2.at(y + 1, x)) + - weight_eps; - if (!mask1.at(y, x) || !mask1.at(y + 1, x) || - !mask2.at(y, x) || !mask2.at(y + 1, x)) - weight += bad_region_penalty_; - graph.addEdges(v, v + img_size.width, weight, weight); - } - } - } -} - - + const Mat &mask1, const Mat &mask2, GCGraph &graph) +{ + const Size img_size = img1.size(); + + // Set terminal weights + for (int y = 0; y < img_size.height; ++y) + { + for (int x = 0; x < img_size.width; ++x) + { + int v = graph.addVtx(); + graph.addTermWeights(v, mask1.at(y, x) ? terminal_cost_ : 0.f, + mask2.at(y, x) ? terminal_cost_ : 0.f); + } + } + + // Set regular edge weights + const float weight_eps = 1.f; + for (int y = 0; y < img_size.height; ++y) + { + for (int x = 0; x < img_size.width; ++x) + { + int v = y * img_size.width + x; + if (x < img_size.width - 1) + { + float weight = normL2(img1.at(y, x), img2.at(y, x)) + + normL2(img1.at(y, x + 1), img2.at(y, x + 1)) + + weight_eps; + if (!mask1.at(y, x) || !mask1.at(y, x + 1) || + !mask2.at(y, x) || !mask2.at(y, x + 1)) + weight += bad_region_penalty_; + graph.addEdges(v, v + 1, weight, weight); + } + if (y < img_size.height - 1) + { + float weight = normL2(img1.at(y, x), img2.at(y, x)) + + normL2(img1.at(y + 1, x), img2.at(y + 1, x)) + + weight_eps; + if (!mask1.at(y, x) || !mask1.at(y + 1, x) || + !mask2.at(y, x) || !mask2.at(y + 1, x)) + weight += bad_region_penalty_; + graph.addEdges(v, v + img_size.width, weight, weight); + } + } + } +} + + void GraphCutSeamFinder::Impl::setGraphWeightsColorGrad( - const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, - const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, - GCGraph &graph) -{ - const Size img_size = img1.size(); - - // Set terminal weights - for (int y = 0; y < img_size.height; ++y) - { - for (int x = 0; x < img_size.width; ++x) - { - int v = graph.addVtx(); - graph.addTermWeights(v, mask1.at(y, x) ? terminal_cost_ : 0.f, - mask2.at(y, x) ? terminal_cost_ : 0.f); - } - } - - // Set regular edge weights - const float weight_eps = 1e-3f; - for (int y = 0; y < img_size.height; ++y) - { - for (int x = 0; x < img_size.width; ++x) - { - int v = y * img_size.width + x; - if (x < img_size.width - 1) - { + const Mat &img1, const Mat &img2, const Mat &dx1, const Mat &dx2, + const Mat &dy1, const Mat &dy2, const Mat &mask1, const Mat &mask2, + GCGraph &graph) +{ + const Size img_size = img1.size(); + + // Set terminal weights + for (int y = 0; y < img_size.height; ++y) + { + for (int x = 0; x < img_size.width; ++x) + { + int v = graph.addVtx(); + graph.addTermWeights(v, mask1.at(y, x) ? terminal_cost_ : 0.f, + mask2.at(y, x) ? terminal_cost_ : 0.f); + } + } + + // Set regular edge weights + const float weight_eps = 1.f; + for (int y = 0; y < img_size.height; ++y) + { + for (int x = 0; x < img_size.width; ++x) + { + int v = y * img_size.width + x; + if (x < img_size.width - 1) + { float grad = dx1.at(y, x) + dx1.at(y, x + 1) + - dx2.at(y, x) + dx2.at(y, x + 1) + weight_eps; - float weight = (normL2(img1.at(y, x), img2.at(y, x)) + - normL2(img1.at(y, x + 1), img2.at(y, x + 1))) / grad + - weight_eps; - if (!mask1.at(y, x) || !mask1.at(y, x + 1) || - !mask2.at(y, x) || !mask2.at(y, x + 1)) - weight += bad_region_penalty_; - graph.addEdges(v, v + 1, weight, weight); - } - if (y < img_size.height - 1) - { + dx2.at(y, x) + dx2.at(y, x + 1) + weight_eps; + float weight = (normL2(img1.at(y, x), img2.at(y, x)) + + normL2(img1.at(y, x + 1), img2.at(y, x + 1))) / grad + + weight_eps; + if (!mask1.at(y, x) || !mask1.at(y, x + 1) || + !mask2.at(y, x) || !mask2.at(y, x + 1)) + weight += bad_region_penalty_; + graph.addEdges(v, v + 1, weight, weight); + } + if (y < img_size.height - 1) + { float grad = dy1.at(y, x) + dy1.at(y + 1, x) + - dy2.at(y, x) + dy2.at(y + 1, x) + weight_eps; - float weight = (normL2(img1.at(y, x), img2.at(y, x)) + - normL2(img1.at(y + 1, x), img2.at(y + 1, x))) / grad + - weight_eps; - if (!mask1.at(y, x) || !mask1.at(y + 1, x) || - !mask2.at(y, x) || !mask2.at(y + 1, x)) - weight += bad_region_penalty_; - graph.addEdges(v, v + img_size.width, weight, weight); - } - } - } -} - - -void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi) -{ - Mat img1 = images_[first], img2 = images_[second]; - Mat dx1 = dx_[first], dx2 = dx_[second]; - Mat dy1 = dy_[first], dy2 = dy_[second]; - Mat mask1 = masks_[first], mask2 = masks_[second]; - Point tl1 = corners_[first], tl2 = corners_[second]; - - const int gap = 10; - Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); - Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); - Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); - Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); - Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); - Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); - Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); - Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); - - // Cut subimages and submasks with some gap - for (int y = -gap; y < roi.height + gap; ++y) - { - for (int x = -gap; x < roi.width + gap; ++x) - { - int y1 = roi.y - tl1.y + y; - int x1 = roi.x - tl1.x + x; - if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) - { - subimg1.at(y + gap, x + gap) = img1.at(y1, x1); - submask1.at(y + gap, x + gap) = mask1.at(y1, x1); - subdx1.at(y + gap, x + gap) = dx1.at(y1, x1); - subdy1.at(y + gap, x + gap) = dy1.at(y1, x1); - } - else - { - subimg1.at(y + gap, x + gap) = Point3f(0, 0, 0); - submask1.at(y + gap, x + gap) = 0; - subdx1.at(y + gap, x + gap) = 0.f; - subdy1.at(y + gap, x + gap) = 0.f; - } - - int y2 = roi.y - tl2.y + y; - int x2 = roi.x - tl2.x + x; - if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) - { - subimg2.at(y + gap, x + gap) = img2.at(y2, x2); - submask2.at(y + gap, x + gap) = mask2.at(y2, x2); - subdx2.at(y + gap, x + gap) = dx2.at(y2, x2); - subdy2.at(y + gap, x + gap) = dy2.at(y2, x2); - } - else - { - subimg2.at(y + gap, x + gap) = Point3f(0, 0, 0); - submask2.at(y + gap, x + gap) = 0; - subdx2.at(y + gap, x + gap) = 0.f; - subdy2.at(y + gap, x + gap) = 0.f; - } - } - } - - const int vertex_count = (roi.height + 2 * gap) * (roi.width + 2 * gap); - const int edge_count = (roi.height - 1 + 2 * gap) * (roi.width + 2 * gap) + - (roi.width - 1 + 2 * gap) * (roi.height + 2 * gap); - GCGraph graph(vertex_count, edge_count); - - switch (cost_type_) - { - case GraphCutSeamFinder::COST_COLOR: - setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph); - break; - case GraphCutSeamFinder::COST_COLOR_GRAD: - setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2, - submask1, submask2, graph); - break; - default: - CV_Error(CV_StsBadArg, "unsupported pixel similarity measure"); - } - - graph.maxFlow(); - - for (int y = 0; y < roi.height; ++y) - { - for (int x = 0; x < roi.width; ++x) - { - if (graph.inSourceSegment((y + gap) * (roi.width + 2 * gap) + x + gap)) - { - if (mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x)) - mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; - } - else - { - if (mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x)) - mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; - } - } - } -} - - -GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty) - : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {} - - -void GraphCutSeamFinder::find(const vector &src, const vector &corners, - vector &masks) -{ - impl_->find(src, corners, masks); -} + dy2.at(y, x) + dy2.at(y + 1, x) + weight_eps; + float weight = (normL2(img1.at(y, x), img2.at(y, x)) + + normL2(img1.at(y + 1, x), img2.at(y + 1, x))) / grad + + weight_eps; + if (!mask1.at(y, x) || !mask1.at(y + 1, x) || + !mask2.at(y, x) || !mask2.at(y + 1, x)) + weight += bad_region_penalty_; + graph.addEdges(v, v + img_size.width, weight, weight); + } + } + } +} + + +void GraphCutSeamFinder::Impl::findInPair(size_t first, size_t second, Rect roi) +{ + Mat img1 = images_[first], img2 = images_[second]; + Mat dx1 = dx_[first], dx2 = dx_[second]; + Mat dy1 = dy_[first], dy2 = dy_[second]; + Mat mask1 = masks_[first], mask2 = masks_[second]; + Point tl1 = corners_[first], tl2 = corners_[second]; + + const int gap = 10; + Mat subimg1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); + Mat subimg2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32FC3); + Mat submask1(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); + Mat submask2(roi.height + 2 * gap, roi.width + 2 * gap, CV_8U); + Mat subdx1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); + Mat subdy1(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); + Mat subdx2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); + Mat subdy2(roi.height + 2 * gap, roi.width + 2 * gap, CV_32F); + + // Cut subimages and submasks with some gap + for (int y = -gap; y < roi.height + gap; ++y) + { + for (int x = -gap; x < roi.width + gap; ++x) + { + int y1 = roi.y - tl1.y + y; + int x1 = roi.x - tl1.x + x; + if (y1 >= 0 && x1 >= 0 && y1 < img1.rows && x1 < img1.cols) + { + subimg1.at(y + gap, x + gap) = img1.at(y1, x1); + submask1.at(y + gap, x + gap) = mask1.at(y1, x1); + subdx1.at(y + gap, x + gap) = dx1.at(y1, x1); + subdy1.at(y + gap, x + gap) = dy1.at(y1, x1); + } + else + { + subimg1.at(y + gap, x + gap) = Point3f(0, 0, 0); + submask1.at(y + gap, x + gap) = 0; + subdx1.at(y + gap, x + gap) = 0.f; + subdy1.at(y + gap, x + gap) = 0.f; + } + + int y2 = roi.y - tl2.y + y; + int x2 = roi.x - tl2.x + x; + if (y2 >= 0 && x2 >= 0 && y2 < img2.rows && x2 < img2.cols) + { + subimg2.at(y + gap, x + gap) = img2.at(y2, x2); + submask2.at(y + gap, x + gap) = mask2.at(y2, x2); + subdx2.at(y + gap, x + gap) = dx2.at(y2, x2); + subdy2.at(y + gap, x + gap) = dy2.at(y2, x2); + } + else + { + subimg2.at(y + gap, x + gap) = Point3f(0, 0, 0); + submask2.at(y + gap, x + gap) = 0; + subdx2.at(y + gap, x + gap) = 0.f; + subdy2.at(y + gap, x + gap) = 0.f; + } + } + } + + const int vertex_count = (roi.height + 2 * gap) * (roi.width + 2 * gap); + const int edge_count = (roi.height - 1 + 2 * gap) * (roi.width + 2 * gap) + + (roi.width - 1 + 2 * gap) * (roi.height + 2 * gap); + GCGraph graph(vertex_count, edge_count); + + switch (cost_type_) + { + case GraphCutSeamFinder::COST_COLOR: + setGraphWeightsColor(subimg1, subimg2, submask1, submask2, graph); + break; + case GraphCutSeamFinder::COST_COLOR_GRAD: + setGraphWeightsColorGrad(subimg1, subimg2, subdx1, subdx2, subdy1, subdy2, + submask1, submask2, graph); + break; + default: + CV_Error(CV_StsBadArg, "unsupported pixel similarity measure"); + } + + graph.maxFlow(); + + for (int y = 0; y < roi.height; ++y) + { + for (int x = 0; x < roi.width; ++x) + { + if (graph.inSourceSegment((y + gap) * (roi.width + 2 * gap) + x + gap)) + { + if (mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x)) + mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x) = 0; + } + else + { + if (mask2.at(roi.y - tl2.y + y, roi.x - tl2.x + x)) + mask1.at(roi.y - tl1.y + y, roi.x - tl1.x + x) = 0; + } + } + } +} + + +GraphCutSeamFinder::GraphCutSeamFinder(int cost_type, float terminal_cost, float bad_region_penalty) + : impl_(new Impl(cost_type, terminal_cost, bad_region_penalty)) {} + + +void GraphCutSeamFinder::find(const vector &src, const vector &corners, + vector &masks) +{ + impl_->find(src, corners, masks); +} diff --git a/modules/stitching/warpers.hpp b/modules/stitching/warpers.hpp index 51d6209..597ee4f 100644 --- a/modules/stitching/warpers.hpp +++ b/modules/stitching/warpers.hpp @@ -38,119 +38,122 @@ // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // -//M*/ -#ifndef __OPENCV_WARPERS_HPP__ -#define __OPENCV_WARPERS_HPP__ - -#include "precomp.hpp" - -class Warper -{ -public: - enum { PLANE, CYLINDRICAL, SPHERICAL }; - static cv::Ptr createByCameraFocal(float focal, int type); - - virtual ~Warper() {} - virtual cv::Point warp(const cv::Mat &src, float focal, const cv::Mat& R, cv::Mat &dst, - int interp_mode = cv::INTER_LINEAR, int border_mode = cv::BORDER_REFLECT) = 0; -}; - - -struct ProjectorBase -{ - void setTransformation(const cv::Mat& R); - - cv::Size size; - float focal; - float r[9]; - float rinv[9]; - float scale; -}; - - -template -class WarperBase : public Warper -{ -public: - cv::Point warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst, - int interp_mode, int border_mode); - -protected: - // Detects ROI of the destination image. It's correct for any projection. - virtual void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); - - // Detects ROI of the destination image by walking over image border. - // Correctness for any projection isn't guaranteed. - void detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br); - - cv::Size src_size_; - P projector_; -}; - - -struct PlaneProjector : ProjectorBase -{ - void mapForward(float x, float y, float &u, float &v); - void mapBackward(float u, float v, float &x, float &y); - - float plane_dist; -}; - - -// Projects image onto z = plane_dist plane -class PlaneWarper : public WarperBase -{ -public: - PlaneWarper(float plane_dist = 1.f, float scale = 1.f) - { - projector_.plane_dist = plane_dist; - projector_.scale = scale; - } - -private: - void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); -}; - - -struct SphericalProjector : ProjectorBase -{ - void mapForward(float x, float y, float &u, float &v); - void mapBackward(float u, float v, float &x, float &y); -}; - - -// Projects image onto unit sphere with origin at (0, 0, 0). -// Poles are located at (0, -1, 0) and (0, 1, 0) points. -class SphericalWarper : public WarperBase -{ -public: - SphericalWarper(float scale = 300.f) { projector_.scale = scale; } - -private: - void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); -}; - - -struct CylindricalProjector : ProjectorBase -{ - void mapForward(float x, float y, float &u, float &v); - void mapBackward(float u, float v, float &x, float &y); -}; - - -// Projects image onto x * x + z * z = 1 cylinder -class CylindricalWarper : public WarperBase -{ -public: - CylindricalWarper(float scale = 300.f) { projector_.scale = scale; } - -private: - void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br) - { - WarperBase::detectResultRoiByBorder(dst_tl, dst_br); - } -}; - -#include "warpers_inl.hpp" - -#endif // __OPENCV_WARPERS_HPP__ +//M*/ +#ifndef __OPENCV_WARPERS_HPP__ +#define __OPENCV_WARPERS_HPP__ + +#include "precomp.hpp" + +class Warper +{ +public: + enum { PLANE, CYLINDRICAL, SPHERICAL }; + static cv::Ptr createByCameraFocal(float focal, int type); + + virtual ~Warper() {} + virtual cv::Point warp(const cv::Mat &src, float focal, const cv::Mat& R, cv::Mat &dst, + int interp_mode = cv::INTER_LINEAR, int border_mode = cv::BORDER_REFLECT) = 0; + virtual cv::Rect warpRoi(const cv::Size &sz, float focal, const cv::Mat &R) = 0; +}; + + +struct ProjectorBase +{ + void setTransformation(const cv::Mat& R); + + cv::Size size; + float focal; + float r[9]; + float rinv[9]; + float scale; +}; + + +template +class WarperBase : public Warper +{ +public: + cv::Point warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst, + int interp_mode, int border_mode); + + cv::Rect warpRoi(const cv::Size &sz, float focal, const cv::Mat &R); + +protected: + // Detects ROI of the destination image. It's correct for any projection. + virtual void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); + + // Detects ROI of the destination image by walking over image border. + // Correctness for any projection isn't guaranteed. + void detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br); + + cv::Size src_size_; + P projector_; +}; + + +struct PlaneProjector : ProjectorBase +{ + void mapForward(float x, float y, float &u, float &v); + void mapBackward(float u, float v, float &x, float &y); + + float plane_dist; +}; + + +// Projects image onto z = plane_dist plane +class PlaneWarper : public WarperBase +{ +public: + PlaneWarper(float plane_dist = 1.f, float scale = 1.f) + { + projector_.plane_dist = plane_dist; + projector_.scale = scale; + } + +private: + void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); +}; + + +struct SphericalProjector : ProjectorBase +{ + void mapForward(float x, float y, float &u, float &v); + void mapBackward(float u, float v, float &x, float &y); +}; + + +// Projects image onto unit sphere with origin at (0, 0, 0). +// Poles are located at (0, -1, 0) and (0, 1, 0) points. +class SphericalWarper : public WarperBase +{ +public: + SphericalWarper(float scale = 300.f) { projector_.scale = scale; } + +private: + void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br); +}; + + +struct CylindricalProjector : ProjectorBase +{ + void mapForward(float x, float y, float &u, float &v); + void mapBackward(float u, float v, float &x, float &y); +}; + + +// Projects image onto x * x + z * z = 1 cylinder +class CylindricalWarper : public WarperBase +{ +public: + CylindricalWarper(float scale = 300.f) { projector_.scale = scale; } + +private: + void detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br) + { + WarperBase::detectResultRoiByBorder(dst_tl, dst_br); + } +}; + +#include "warpers_inl.hpp" + +#endif // __OPENCV_WARPERS_HPP__ diff --git a/modules/stitching/warpers_inl.hpp b/modules/stitching/warpers_inl.hpp index d55558c..e385ff1 100644 --- a/modules/stitching/warpers_inl.hpp +++ b/modules/stitching/warpers_inl.hpp @@ -38,202 +38,218 @@ // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // -//M*/ -#ifndef __OPENCV_WARPERS_INL_HPP__ -#define __OPENCV_WARPERS_INL_HPP__ - -#include "warpers.hpp" // Make your IDE see declarations - -template -cv::Point WarperBase

::warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst, - int interp_mode, int border_mode) -{ - src_size_ = src.size(); - - projector_.size = src.size(); - projector_.focal = focal; - projector_.setTransformation(R); - - cv::Point dst_tl, dst_br; - detectResultRoi(dst_tl, dst_br); - - cv::Mat xmap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F); - cv::Mat ymap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F); - - float x, y; - for (int v = dst_tl.y; v <= dst_br.y; ++v) - { - for (int u = dst_tl.x; u <= dst_br.x; ++u) - { - projector_.mapBackward(static_cast(u), static_cast(v), x, y); - xmap.at(v - dst_tl.y, u - dst_tl.x) = x; - ymap.at(v - dst_tl.y, u - dst_tl.x) = y; - } - } - - dst.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type()); - remap(src, dst, xmap, ymap, interp_mode, border_mode); - - return dst_tl; -} - - -template -void WarperBase

::detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br) -{ - float tl_uf = std::numeric_limits::max(); - float tl_vf = std::numeric_limits::max(); - float br_uf = -std::numeric_limits::max(); - float br_vf = -std::numeric_limits::max(); - - float u, v; - for (int y = 0; y < src_size_.height; ++y) - { - for (int x = 0; x < src_size_.width; ++x) - { - projector_.mapForward(static_cast(x), static_cast(y), u, v); - tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); - br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); - } - } - - dst_tl.x = static_cast(tl_uf); - dst_tl.y = static_cast(tl_vf); - dst_br.x = static_cast(br_uf); - dst_br.y = static_cast(br_vf); -} - - -template -void WarperBase

::detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br) -{ - float tl_uf = std::numeric_limits::max(); - float tl_vf = std::numeric_limits::max(); - float br_uf = -std::numeric_limits::max(); - float br_vf = -std::numeric_limits::max(); - - float u, v; - for (float x = 0; x < src_size_.width; ++x) - { - projector_.mapForward(static_cast(x), 0, u, v); - tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); - br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); - - projector_.mapForward(static_cast(x), static_cast(src_size_.height - 1), u, v); - tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); - br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); - } - for (int y = 0; y < src_size_.height; ++y) - { - projector_.mapForward(0, static_cast(y), u, v); - tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); - br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); - - projector_.mapForward(static_cast(src_size_.width - 1), static_cast(y), u, v); - tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); - br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); - } - - dst_tl.x = static_cast(tl_uf); - dst_tl.y = static_cast(tl_vf); - dst_br.x = static_cast(br_uf); - dst_br.y = static_cast(br_vf); -} - - -inline -void PlaneProjector::mapForward(float x, float y, float &u, float &v) -{ - x -= size.width * 0.5f; - y -= size.height * 0.5f; - - float x_ = r[0] * x + r[1] * y + r[2] * focal; - float y_ = r[3] * x + r[4] * y + r[5] * focal; - float z_ = r[6] * x + r[7] * y + r[8] * focal; - - u = scale * x_ / z_ * plane_dist; - v = scale * y_ / z_ * plane_dist; -} - - -inline -void PlaneProjector::mapBackward(float u, float v, float &x, float &y) -{ - float x_ = u / scale; - float y_ = v / scale; - - float z; - x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * plane_dist; - y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * plane_dist; - z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * plane_dist; - - x = focal * x / z + size.width * 0.5f; - y = focal * y / z + size.height * 0.5f; -} - - -inline -void SphericalProjector::mapForward(float x, float y, float &u, float &v) -{ - x -= size.width * 0.5f; - y -= size.height * 0.5f; - - float x_ = r[0] * x + r[1] * y + r[2] * focal; - float y_ = r[3] * x + r[4] * y + r[5] * focal; - float z_ = r[6] * x + r[7] * y + r[8] * focal; - - u = scale * atan2f(x_, z_); - v = scale * (static_cast(CV_PI) - acosf(y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_))); -} - - -inline -void SphericalProjector::mapBackward(float u, float v, float &x, float &y) -{ - float sinv = sinf(static_cast(CV_PI) - v / scale); - float x_ = sinv * sinf(u / scale); - float y_ = cosf(static_cast(CV_PI) - v / scale); - float z_ = sinv * cosf(u / scale); - - float z; - x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_; - y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_; - z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_; - - x = focal * x / z + size.width * 0.5f; - y = focal * y / z + size.height * 0.5f; -} - - -inline -void CylindricalProjector::mapForward(float x, float y, float &u, float &v) -{ - x -= size.width * 0.5f; - y -= size.height * 0.5f; - - float x_ = r[0] * x + r[1] * y + r[2] * focal; - float y_ = r[3] * x + r[4] * y + r[5] * focal; - float z_ = r[6] * x + r[7] * y + r[8] * focal; - - u = scale * atan2f(x_, z_); - v = scale * y_ / sqrtf(x_ * x_ + z_ * z_); -} - - -inline -void CylindricalProjector::mapBackward(float u, float v, float &x, float &y) -{ - float x_ = sinf(u / scale); - float y_ = v / scale; - float z_ = cosf(u / scale); - - float z; - x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_; - y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_; - z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_; - - x = focal * x / z + size.width * 0.5f; - y = focal * y / z + size.height * 0.5f; -} - -#endif // __OPENCV_WARPERS_INL_HPP__ +//M*/ +#ifndef __OPENCV_WARPERS_INL_HPP__ +#define __OPENCV_WARPERS_INL_HPP__ + +#include "warpers.hpp" // Make your IDE see declarations + +template +cv::Point WarperBase

::warp(const cv::Mat &src, float focal, const cv::Mat &R, cv::Mat &dst, + int interp_mode, int border_mode) +{ + src_size_ = src.size(); + + projector_.size = src.size(); + projector_.focal = focal; + projector_.setTransformation(R); + + cv::Point dst_tl, dst_br; + detectResultRoi(dst_tl, dst_br); + + cv::Mat xmap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F); + cv::Mat ymap(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, CV_32F); + + float x, y; + for (int v = dst_tl.y; v <= dst_br.y; ++v) + { + for (int u = dst_tl.x; u <= dst_br.x; ++u) + { + projector_.mapBackward(static_cast(u), static_cast(v), x, y); + xmap.at(v - dst_tl.y, u - dst_tl.x) = x; + ymap.at(v - dst_tl.y, u - dst_tl.x) = y; + } + } + + dst.create(dst_br.y - dst_tl.y + 1, dst_br.x - dst_tl.x + 1, src.type()); + remap(src, dst, xmap, ymap, interp_mode, border_mode); + + return dst_tl; +} + + +template +cv::Rect WarperBase

::warpRoi(const cv::Size &sz, float focal, const cv::Mat &R) +{ + src_size_ = sz; + + projector_.size = sz; + projector_.focal = focal; + projector_.setTransformation(R); + + cv::Point dst_tl, dst_br; + detectResultRoi(dst_tl, dst_br); + + return cv::Rect(dst_tl, cv::Point(dst_br.x + 1, dst_br.y + 1)); +} + + +template +void WarperBase

::detectResultRoi(cv::Point &dst_tl, cv::Point &dst_br) +{ + float tl_uf = std::numeric_limits::max(); + float tl_vf = std::numeric_limits::max(); + float br_uf = -std::numeric_limits::max(); + float br_vf = -std::numeric_limits::max(); + + float u, v; + for (int y = 0; y < src_size_.height; ++y) + { + for (int x = 0; x < src_size_.width; ++x) + { + projector_.mapForward(static_cast(x), static_cast(y), u, v); + tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); + br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); + } + } + + dst_tl.x = static_cast(tl_uf); + dst_tl.y = static_cast(tl_vf); + dst_br.x = static_cast(br_uf); + dst_br.y = static_cast(br_vf); +} + + +template +void WarperBase

::detectResultRoiByBorder(cv::Point &dst_tl, cv::Point &dst_br) +{ + float tl_uf = std::numeric_limits::max(); + float tl_vf = std::numeric_limits::max(); + float br_uf = -std::numeric_limits::max(); + float br_vf = -std::numeric_limits::max(); + + float u, v; + for (float x = 0; x < src_size_.width; ++x) + { + projector_.mapForward(static_cast(x), 0, u, v); + tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); + br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); + + projector_.mapForward(static_cast(x), static_cast(src_size_.height - 1), u, v); + tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); + br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); + } + for (int y = 0; y < src_size_.height; ++y) + { + projector_.mapForward(0, static_cast(y), u, v); + tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); + br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); + + projector_.mapForward(static_cast(src_size_.width - 1), static_cast(y), u, v); + tl_uf = std::min(tl_uf, u); tl_vf = std::min(tl_vf, v); + br_uf = std::max(br_uf, u); br_vf = std::max(br_vf, v); + } + + dst_tl.x = static_cast(tl_uf); + dst_tl.y = static_cast(tl_vf); + dst_br.x = static_cast(br_uf); + dst_br.y = static_cast(br_vf); +} + + +inline +void PlaneProjector::mapForward(float x, float y, float &u, float &v) +{ + x -= size.width * 0.5f; + y -= size.height * 0.5f; + + float x_ = r[0] * x + r[1] * y + r[2] * focal; + float y_ = r[3] * x + r[4] * y + r[5] * focal; + float z_ = r[6] * x + r[7] * y + r[8] * focal; + + u = scale * x_ / z_ * plane_dist; + v = scale * y_ / z_ * plane_dist; +} + + +inline +void PlaneProjector::mapBackward(float u, float v, float &x, float &y) +{ + float x_ = u / scale; + float y_ = v / scale; + + float z; + x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * plane_dist; + y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * plane_dist; + z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * plane_dist; + + x = focal * x / z + size.width * 0.5f; + y = focal * y / z + size.height * 0.5f; +} + + +inline +void SphericalProjector::mapForward(float x, float y, float &u, float &v) +{ + x -= size.width * 0.5f; + y -= size.height * 0.5f; + + float x_ = r[0] * x + r[1] * y + r[2] * focal; + float y_ = r[3] * x + r[4] * y + r[5] * focal; + float z_ = r[6] * x + r[7] * y + r[8] * focal; + + u = scale * atan2f(x_, z_); + v = scale * (static_cast(CV_PI) - acosf(y_ / sqrtf(x_ * x_ + y_ * y_ + z_ * z_))); +} + + +inline +void SphericalProjector::mapBackward(float u, float v, float &x, float &y) +{ + float sinv = sinf(static_cast(CV_PI) - v / scale); + float x_ = sinv * sinf(u / scale); + float y_ = cosf(static_cast(CV_PI) - v / scale); + float z_ = sinv * cosf(u / scale); + + float z; + x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_; + y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_; + z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_; + + x = focal * x / z + size.width * 0.5f; + y = focal * y / z + size.height * 0.5f; +} + + +inline +void CylindricalProjector::mapForward(float x, float y, float &u, float &v) +{ + x -= size.width * 0.5f; + y -= size.height * 0.5f; + + float x_ = r[0] * x + r[1] * y + r[2] * focal; + float y_ = r[3] * x + r[4] * y + r[5] * focal; + float z_ = r[6] * x + r[7] * y + r[8] * focal; + + u = scale * atan2f(x_, z_); + v = scale * y_ / sqrtf(x_ * x_ + z_ * z_); +} + + +inline +void CylindricalProjector::mapBackward(float u, float v, float &x, float &y) +{ + float x_ = sinf(u / scale); + float y_ = v / scale; + float z_ = cosf(u / scale); + + float z; + x = rinv[0] * x_ + rinv[1] * y_ + rinv[2] * z_; + y = rinv[3] * x_ + rinv[4] * y_ + rinv[5] * z_; + z = rinv[6] * x_ + rinv[7] * y_ + rinv[8] * z_; + + x = focal * x / z + size.width * 0.5f; + y = focal * y / z + size.height * 0.5f; +} + +#endif // __OPENCV_WARPERS_INL_HPP__ diff --git a/samples/gpu/driver_api_stereo_multi.cpp b/samples/gpu/driver_api_stereo_multi.cpp index 9f65c46..ce2925b 100644 --- a/samples/gpu/driver_api_stereo_multi.cpp +++ b/samples/gpu/driver_api_stereo_multi.cpp @@ -32,7 +32,11 @@ int main() #include #include +#include +#include #include "opencv2/core/internal.hpp" // For TBB wrappers +#include "tbb/tbb.h" +#include "tbb/mutex.h" using namespace std; using namespace cv; @@ -54,7 +58,7 @@ inline void safeCall_(int code, const char* expr, const char* file, int line) } // Each GPU is associated with its own context -CUcontext contexts[2]; +CUcontext contexts[/*2*/1]; void inline contextOn(int id) { @@ -76,6 +80,10 @@ GpuMat d_result[2]; // CPU result Mat result; +int some[2]; + +tbb::mutex mutex; + int main(int argc, char** argv) { if (argc < 3) @@ -85,11 +93,11 @@ int main(int argc, char** argv) } int num_devices = getCudaEnabledDeviceCount(); - if (num_devices < 2) - { - std::cout << "Two or more GPUs are required\n"; - return -1; - } +// if (num_devices < 2) +// { +// std::cout << "Two or more GPUs are required\n"; +// return -1; +// } for (int i = 0; i < num_devices; ++i) { @@ -123,13 +131,14 @@ int main(int argc, char** argv) // Create context for GPU #0 CUdevice device; safeCall(cuDeviceGet(&device, 0)); - safeCall(cuCtxCreate(&contexts[0], 0, device)); + safeCall(cuGLCtxCreate(&contexts[0], 0, device)); + //safeCall(cuCtxCreate(&contexts[0], 0, device)); contextOff(); - // Create context for GPU #1 - safeCall(cuDeviceGet(&device, 1)); - safeCall(cuCtxCreate(&contexts[1], 0, device)); - contextOff(); +// // Create context for GPU #1 +// safeCall(cuDeviceGet(&device, 0)); +// safeCall(cuCtxCreate(&contexts[1], 0, device)); +// contextOff(); // Split source images for processing on GPU #0 contextOn(0); @@ -139,15 +148,20 @@ int main(int argc, char** argv) contextOff(); // Split source images for processing on the GPU #1 - contextOn(1); + contextOn(0); d_left[1].upload(left.rowRange(left.rows / 2, left.rows)); d_right[1].upload(right.rowRange(right.rows / 2, right.rows)); bm[1] = new StereoBM_GPU(); contextOff(); + some[0] = some[1] = 0; // Execute calculation in two threads using two GPUs - int devices[] = {0, 1}; - parallel_do(devices, devices + 2, Worker()); + vector devices; + for (int i = 0; i < 4; ++i) + devices.push_back(rand()%2); + tbb::parallel_do(&devices[0], &devices[devices.size() - 1], Worker()); + + cout << some[0] << " " << some[1] << endl; // Release the first GPU resources contextOn(0); @@ -159,7 +173,7 @@ int main(int argc, char** argv) contextOff(); // Release the second GPU resources - contextOn(1); + contextOn(0); imshow("GPU #1 result", Mat(d_result[1])); d_left[1].release(); d_right[1].release(); @@ -175,7 +189,9 @@ int main(int argc, char** argv) void Worker::operator()(int device_id) const { - contextOn(device_id); + mutex.lock(); + + contextOn(0); bm[device_id]->operator()(d_left[device_id], d_right[device_id], d_result[device_id]); @@ -184,13 +200,16 @@ void Worker::operator()(int device_id) const << "): finished\n"; contextOff(); + + mutex.unlock(); } void destroyContexts() { safeCall(cuCtxDestroy(contexts[0])); - safeCall(cuCtxDestroy(contexts[1])); + //safeCall(cuCtxDestroy(contexts[1])); } #endif + -- 2.7.4