stitcher.cpp 15.5 KB
Newer Older
wester committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                          License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// 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 "precomp.hpp"

wester committed
45 46
using namespace std;

wester committed
47 48 49 50 51 52 53 54 55 56 57
namespace cv {

Stitcher Stitcher::createDefault(bool try_use_gpu)
{
    Stitcher stitcher;
    stitcher.setRegistrationResol(0.6);
    stitcher.setSeamEstimationResol(0.1);
    stitcher.setCompositingResol(ORIG_RESOL);
    stitcher.setPanoConfidenceThresh(1);
    stitcher.setWaveCorrection(true);
    stitcher.setWaveCorrectKind(detail::WAVE_CORRECT_HORIZ);
wester committed
58 59
    stitcher.setFeaturesMatcher(new detail::BestOf2NearestMatcher(try_use_gpu));
    stitcher.setBundleAdjuster(new detail::BundleAdjusterRay());
wester committed
60

wester committed
61 62
#if defined(HAVE_OPENCV_GPU) && !defined(DYNAMIC_CUDA_SUPPORT)
    if (try_use_gpu && gpu::getCudaEnabledDeviceCount() > 0)
wester committed
63
    {
wester committed
64 65
#if defined(HAVE_OPENCV_NONFREE)
        stitcher.setFeaturesFinder(new detail::SurfFeaturesFinderGpu());
wester committed
66
#else
wester committed
67
        stitcher.setFeaturesFinder(new detail::OrbFeaturesFinder());
wester committed
68
#endif
wester committed
69 70
        stitcher.setWarper(new SphericalWarperGpu());
        stitcher.setSeamFinder(new detail::GraphCutSeamFinderGpu());
wester committed
71 72 73 74
    }
    else
#endif
    {
wester committed
75 76
#ifdef HAVE_OPENCV_NONFREE
        stitcher.setFeaturesFinder(new detail::SurfFeaturesFinder());
wester committed
77
#else
wester committed
78
        stitcher.setFeaturesFinder(new detail::OrbFeaturesFinder());
wester committed
79
#endif
wester committed
80 81
        stitcher.setWarper(new SphericalWarper());
        stitcher.setSeamFinder(new detail::GraphCutSeamFinder(detail::GraphCutSeamFinderBase::COST_COLOR));
wester committed
82 83
    }

wester committed
84 85
    stitcher.setExposureCompensator(new detail::BlocksGainCompensator());
    stitcher.setBlender(new detail::MultiBandBlender(try_use_gpu));
wester committed
86 87 88 89 90

    return stitcher;
}


wester committed
91
Stitcher::Status Stitcher::estimateTransform(InputArray images)
wester committed
92
{
wester committed
93
    return estimateTransform(images, vector<vector<Rect> >());
wester committed
94 95 96
}


wester committed
97
Stitcher::Status Stitcher::estimateTransform(InputArray images, const vector<vector<Rect> > &rois)
wester committed
98
{
wester committed
99
    images.getMatVector(imgs_);
wester committed
100 101 102 103 104 105 106
    rois_ = rois;

    Status status;

    if ((status = matchImages()) != OK)
        return status;

wester committed
107
    estimateCameraParams();
wester committed
108 109 110 111 112 113 114 115

    return OK;
}



Stitcher::Status Stitcher::composePanorama(OutputArray pano)
{
wester committed
116
    return composePanorama(vector<Mat>(), pano);
wester committed
117 118 119
}


wester committed
120
Stitcher::Status Stitcher::composePanorama(InputArray images, OutputArray pano)
wester committed
121 122 123
{
    LOGLN("Warping images (auxiliary)... ");

wester committed
124 125
    vector<Mat> imgs;
    images.getMatVector(imgs);
wester committed
126 127 128 129
    if (!imgs.empty())
    {
        CV_Assert(imgs.size() == imgs_.size());

wester committed
130
        Mat img;
wester committed
131 132 133 134 135 136 137 138 139
        seam_est_imgs_.resize(imgs.size());

        for (size_t i = 0; i < imgs.size(); ++i)
        {
            imgs_[i] = imgs[i];
            resize(imgs[i], img, Size(), seam_scale_, seam_scale_);
            seam_est_imgs_[i] = img.clone();
        }

wester committed
140 141
        vector<Mat> seam_est_imgs_subset;
        vector<Mat> imgs_subset;
wester committed
142 143 144 145 146 147 148 149 150 151 152

        for (size_t i = 0; i < indices_.size(); ++i)
        {
            imgs_subset.push_back(imgs_[indices_[i]]);
            seam_est_imgs_subset.push_back(seam_est_imgs_[indices_[i]]);
        }

        seam_est_imgs_ = seam_est_imgs_subset;
        imgs_ = imgs_subset;
    }

wester committed
153
    Mat &pano_ = pano.getMatRef();
wester committed
154 155 156 157 158

#if ENABLE_LOG
    int64 t = getTickCount();
#endif

wester committed
159 160 161 162 163
    vector<Point> corners(imgs_.size());
    vector<Mat> masks_warped(imgs_.size());
    vector<Mat> images_warped(imgs_.size());
    vector<Size> sizes(imgs_.size());
    vector<Mat> masks(imgs_.size());
wester committed
164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182

    // Prepare image masks
    for (size_t i = 0; i < imgs_.size(); ++i)
    {
        masks[i].create(seam_est_imgs_[i].size(), CV_8U);
        masks[i].setTo(Scalar::all(255));
    }

    // Warp images and their masks
    Ptr<detail::RotationWarper> w = warper_->create(float(warped_image_scale_ * seam_work_aspect_));
    for (size_t i = 0; i < imgs_.size(); ++i)
    {
        Mat_<float> K;
        cameras_[i].K().convertTo(K, CV_32F);
        K(0,0) *= (float)seam_work_aspect_;
        K(0,2) *= (float)seam_work_aspect_;
        K(1,1) *= (float)seam_work_aspect_;
        K(1,2) *= (float)seam_work_aspect_;

wester committed
183
        corners[i] = w->warp(seam_est_imgs_[i], K, cameras_[i].R, INTER_LINEAR, BORDER_REFLECT, images_warped[i]);
wester committed
184 185 186 187 188
        sizes[i] = images_warped[i].size();

        w->warp(masks[i], K, cameras_[i].R, INTER_NEAREST, BORDER_CONSTANT, masks_warped[i]);
    }

wester committed
189
    vector<Mat> images_warped_f(imgs_.size());
a  
Kai Westerkamp committed
190 191
    for (size_t i = 0; i < imgs_.size(); ++i)
        images_warped[i].convertTo(images_warped_f[i], CV_32F);
wester committed
192 193 194 195

    LOGLN("Warping images, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");

    // Find seams
a  
Kai Westerkamp committed
196
    exposure_comp_->feed(corners, images_warped, masks_warped);
wester committed
197 198 199 200 201 202 203 204 205 206 207 208 209
    seam_finder_->find(images_warped_f, corners, masks_warped);

    // Release unused memory
    seam_est_imgs_.clear();
    images_warped.clear();
    images_warped_f.clear();
    masks.clear();

    LOGLN("Compositing...");
#if ENABLE_LOG
    t = getTickCount();
#endif

wester committed
210 211
    Mat img_warped, img_warped_s;
    Mat dilated_mask, seam_mask, mask, mask_warped;
wester committed
212 213 214 215 216 217 218 219

    //double compose_seam_aspect = 1;
    double compose_work_aspect = 1;
    bool is_blender_prepared = false;

    double compose_scale = 1;
    bool is_compose_scale_set = false;

wester committed
220
    Mat full_img, img;
wester committed
221 222 223 224 225 226 227 228 229
    for (size_t img_idx = 0; img_idx < imgs_.size(); ++img_idx)
    {
        LOGLN("Compositing image #" << indices_[img_idx] + 1);

        // Read image and resize it if necessary
        full_img = imgs_[img_idx];
        if (!is_compose_scale_set)
        {
            if (compose_resol_ > 0)
wester committed
230
                compose_scale = min(1.0, sqrt(compose_resol_ * 1e6 / full_img.size().area()));
wester committed
231 232 233 234 235 236 237
            is_compose_scale_set = true;

            // Compute relative scales
            //compose_seam_aspect = compose_scale / seam_scale_;
            compose_work_aspect = compose_scale / work_scale_;

            // Update warped image scale
a  
Kai Westerkamp committed
238 239
            warped_image_scale_ *= static_cast<float>(compose_work_aspect);
            w = warper_->create((float)warped_image_scale_);
wester committed
240 241 242 243 244

            // Update corners and sizes
            for (size_t i = 0; i < imgs_.size(); ++i)
            {
                // Update intrinsics
a  
Kai Westerkamp committed
245 246 247
                cameras_[i].focal *= compose_work_aspect;
                cameras_[i].ppx *= compose_work_aspect;
                cameras_[i].ppy *= compose_work_aspect;
wester committed
248 249 250 251 252 253 254 255 256 257

                // Update corner and size
                Size sz = full_img_sizes_[i];
                if (std::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);
                }

                Mat K;
a  
Kai Westerkamp committed
258 259
                cameras_[i].K().convertTo(K, CV_32F);
                Rect roi = w->warpRoi(sz, K, cameras_[i].R);
wester committed
260 261 262 263 264 265 266 267 268 269 270 271
                corners[i] = roi.tl();
                sizes[i] = roi.size();
            }
        }
        if (std::abs(compose_scale - 1) > 1e-1)
            resize(full_img, img, Size(), compose_scale, compose_scale);
        else
            img = full_img;
        full_img.release();
        Size img_size = img.size();

        Mat K;
a  
Kai Westerkamp committed
272
        cameras_[img_idx].K().convertTo(K, CV_32F);
wester committed
273 274

        // Warp the current image
wester committed
275
        w->warp(img, K, cameras_[img_idx].R, INTER_LINEAR, BORDER_REFLECT, img_warped);
wester committed
276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293

        // Warp the current image mask
        mask.create(img_size, CV_8U);
        mask.setTo(Scalar::all(255));
        w->warp(mask, K, cameras_[img_idx].R, INTER_NEAREST, BORDER_CONSTANT, mask_warped);

        // Compensate exposure
        exposure_comp_->apply((int)img_idx, corners[img_idx], img_warped, mask_warped);

        img_warped.convertTo(img_warped_s, CV_16S);
        img_warped.release();
        img.release();
        mask.release();

        // Make sure seam mask has proper size
        dilate(masks_warped[img_idx], dilated_mask, Mat());
        resize(dilated_mask, seam_mask, mask_warped.size());

wester committed
294
        mask_warped = seam_mask & mask_warped;
wester committed
295 296 297 298 299 300 301 302 303 304 305

        if (!is_blender_prepared)
        {
            blender_->prepare(corners, sizes);
            is_blender_prepared = true;
        }

        // Blend the current image
        blender_->feed(img_warped_s, mask_warped, corners[img_idx]);
    }

wester committed
306
    Mat result, result_mask;
wester committed
307 308 309 310 311 312
    blender_->blend(result, result_mask);

    LOGLN("Compositing, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");

    // Preliminary result is in CV_16SC3 format, but all values are in [0,255] range,
    // so convert it to avoid user confusing
wester committed
313
    result.convertTo(pano_, CV_8U);
wester committed
314 315 316 317 318

    return OK;
}


wester committed
319
Stitcher::Status Stitcher::stitch(InputArray images, OutputArray pano)
wester committed
320 321 322 323 324 325 326 327
{
    Status status = estimateTransform(images);
    if (status != OK)
        return status;
    return composePanorama(pano);
}


wester committed
328
Stitcher::Status Stitcher::stitch(InputArray images, const vector<vector<Rect> > &rois, OutputArray pano)
wester committed
329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
{
    Status status = estimateTransform(images, rois);
    if (status != OK)
        return status;
    return composePanorama(pano);
}


Stitcher::Status Stitcher::matchImages()
{
    if ((int)imgs_.size() < 2)
    {
        LOGLN("Need more images");
        return ERR_NEED_MORE_IMGS;
    }

    work_scale_ = 1;
    seam_work_aspect_ = 1;
    seam_scale_ = 1;
    bool is_work_scale_set = false;
    bool is_seam_scale_set = false;
wester committed
350
    Mat full_img, img;
wester committed
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
    features_.resize(imgs_.size());
    seam_est_imgs_.resize(imgs_.size());
    full_img_sizes_.resize(imgs_.size());

    LOGLN("Finding features...");
#if ENABLE_LOG
    int64 t = getTickCount();
#endif

    for (size_t i = 0; i < imgs_.size(); ++i)
    {
        full_img = imgs_[i];
        full_img_sizes_[i] = full_img.size();

        if (registr_resol_ < 0)
        {
            img = full_img;
            work_scale_ = 1;
            is_work_scale_set = true;
        }
        else
        {
            if (!is_work_scale_set)
            {
wester committed
375
                work_scale_ = min(1.0, sqrt(registr_resol_ * 1e6 / full_img.size().area()));
wester committed
376 377 378 379 380 381
                is_work_scale_set = true;
            }
            resize(full_img, img, Size(), work_scale_, work_scale_);
        }
        if (!is_seam_scale_set)
        {
wester committed
382
            seam_scale_ = min(1.0, sqrt(seam_est_resol_ * 1e6 / full_img.size().area()));
wester committed
383 384 385 386 387
            seam_work_aspect_ = seam_scale_ / work_scale_;
            is_seam_scale_set = true;
        }

        if (rois_.empty())
a  
Kai Westerkamp committed
388
            (*features_finder_)(img, features_[i]);
wester committed
389 390
        else
        {
wester committed
391
            vector<Rect> rois(rois_[i].size());
wester committed
392 393 394 395
            for (size_t j = 0; j < rois_[i].size(); ++j)
            {
                Point tl(cvRound(rois_[i][j].x * work_scale_), cvRound(rois_[i][j].y * work_scale_));
                Point br(cvRound(rois_[i][j].br().x * work_scale_), cvRound(rois_[i][j].br().y * work_scale_));
a  
Kai Westerkamp committed
396
                rois[j] = Rect(tl, br);
wester committed
397
            }
a  
Kai Westerkamp committed
398
            (*features_finder_)(img, features_[i], rois);
wester committed
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423
        }
        features_[i].img_idx = (int)i;
        LOGLN("Features in image #" << i+1 << ": " << features_[i].keypoints.size());

        resize(full_img, img, Size(), seam_scale_, seam_scale_);
        seam_est_imgs_[i] = img.clone();
    }

    // Do it to save memory
    features_finder_->collectGarbage();
    full_img.release();
    img.release();

    LOGLN("Finding features, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");

    LOG("Pairwise matching");
#if ENABLE_LOG
    t = getTickCount();
#endif
    (*features_matcher_)(features_, pairwise_matches_, matching_mask_);
    features_matcher_->collectGarbage();
    LOGLN("Pairwise matching, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");

    // Leave only images we are sure are from the same panorama
    indices_ = detail::leaveBiggestComponent(features_, pairwise_matches_, (float)conf_thresh_);
wester committed
424 425 426
    vector<Mat> seam_est_imgs_subset;
    vector<Mat> imgs_subset;
    vector<Size> full_img_sizes_subset;
wester committed
427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446
    for (size_t i = 0; i < indices_.size(); ++i)
    {
        imgs_subset.push_back(imgs_[indices_[i]]);
        seam_est_imgs_subset.push_back(seam_est_imgs_[indices_[i]]);
        full_img_sizes_subset.push_back(full_img_sizes_[indices_[i]]);
    }
    seam_est_imgs_ = seam_est_imgs_subset;
    imgs_ = imgs_subset;
    full_img_sizes_ = full_img_sizes_subset;

    if ((int)imgs_.size() < 2)
    {
        LOGLN("Need more images");
        return ERR_NEED_MORE_IMGS;
    }

    return OK;
}


wester committed
447
void Stitcher::estimateCameraParams()
wester committed
448
{
a  
Kai Westerkamp committed
449
    detail::HomographyBasedEstimator estimator;
wester committed
450
    estimator(features_, pairwise_matches_, cameras_);
wester committed
451 452 453 454 455 456

    for (size_t i = 0; i < cameras_.size(); ++i)
    {
        Mat R;
        cameras_[i].R.convertTo(R, CV_32F);
        cameras_[i].R = R;
wester committed
457
        LOGLN("Initial intrinsic parameters #" << indices_[i] + 1 << ":\n " << cameras_[i].K());
wester committed
458 459 460
    }

    bundle_adjuster_->setConfThresh(conf_thresh_);
wester committed
461
    (*bundle_adjuster_)(features_, pairwise_matches_, cameras_);
wester committed
462 463

    // Find median focal length and use it as final image scale
wester committed
464
    vector<double> focals;
wester committed
465 466
    for (size_t i = 0; i < cameras_.size(); ++i)
    {
wester committed
467
        LOGLN("Camera #" << indices_[i] + 1 << ":\n" << cameras_[i].K());
wester committed
468 469 470 471 472 473 474 475 476 477 478
        focals.push_back(cameras_[i].focal);
    }

    std::sort(focals.begin(), focals.end());
    if (focals.size() % 2 == 1)
        warped_image_scale_ = static_cast<float>(focals[focals.size() / 2]);
    else
        warped_image_scale_ = static_cast<float>(focals[focals.size() / 2 - 1] + focals[focals.size() / 2]) * 0.5f;

    if (do_wave_correct_)
    {
wester committed
479
        vector<Mat> rmats;
wester committed
480 481 482 483 484 485 486 487 488
        for (size_t i = 0; i < cameras_.size(); ++i)
            rmats.push_back(cameras_[i].R.clone());
        detail::waveCorrect(rmats, wave_correct_kind_);
        for (size_t i = 0; i < cameras_.size(); ++i)
            cameras_[i].R = rmats[i];
    }
}

} // namespace cv