PageRenderTime 5ms CodeModel.GetById 161ms app.highlight 412ms RepoModel.GetById 181ms app.codeStats 1ms

/toolkits/computer_vision/stitching_detailed.cpp

https://github.com/michaelkook/GraphLab-2
C++ | 737 lines | 622 code | 73 blank | 42 comment | 231 complexity | 15d0d585f262d288b9c8d7f7cfbdc4c1 MD5 | raw file
  1/*  
  2 * Copyright (c) 2009 Carnegie Mellon University. 
  3 *     All rights reserved.
  4 *
  5 *  Licensed under the Apache License, Version 2.0 (the "License");
  6 *  you may not use this file except in compliance with the License.
  7 *  You may obtain a copy of the License at
  8 *
  9 *      http://www.apache.org/licenses/LICENSE-2.0
 10 *
 11 *  Unless required by applicable law or agreed to in writing,
 12 *  software distributed under the License is distributed on an "AS
 13 *  IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
 14 *  express or implied.  See the License for the specific language
 15 *  governing permissions and limitations under the License.
 16 *
 17 * For more about this software visit:
 18 *
 19 *      http://www.graphlab.ml.cmu.edu
 20 *
 21 */
 22
 23
 24
 25#include <fstream>
 26#include <string>
 27#include "opencv2/opencv_modules.hpp"
 28#include "opencv2/highgui/highgui.hpp"
 29#include "opencv2/stitching/detail/autocalib.hpp"
 30#include "opencv2/stitching/detail/blenders.hpp"
 31#include "opencv2/stitching/detail/camera.hpp"
 32#include "opencv2/stitching/detail/exposure_compensate.hpp"
 33#include "opencv2/stitching/detail/matchers.hpp"
 34#include "opencv2/stitching/detail/motion_estimators.hpp"
 35#include "opencv2/stitching/detail/seam_finders.hpp"
 36#include "opencv2/stitching/detail/util.hpp"
 37#include "opencv2/stitching/detail/warpers.hpp"
 38#include "opencv2/stitching/warpers.hpp"
 39
 40using namespace std;
 41using namespace cv;
 42using namespace cv::detail;
 43
 44void printUsage()
 45{
 46    cout <<
 47        "Rotation model images stitcher.\n\n"
 48        "stitching_detailed img1 img2 [...imgN] [flags]\n\n"
 49        "Flags:\n"
 50        "  --preview\n"
 51        "      Run stitching in the preview mode. Works faster than usual mode,\n"
 52        "      but output image will have lower resolution.\n"
 53        "  --try_gpu (yes|no)\n"
 54        "      Try to use GPU. The default value is 'no'. All default values\n"
 55        "      are for CPU mode.\n"
 56        "\nMotion Estimation Flags:\n"
 57        "  --work_megapix <float>\n"
 58        "      Resolution for image registration step. The default is 0.6 Mpx.\n"
 59        "  --features (surf|orb)\n"
 60        "      Type of features used for images matching. The default is surf.\n"
 61        "  --match_conf <float>\n"
 62        "      Confidence for feature matching step. The default is 0.65 for surf and 0.3 for orb.\n"
 63        "  --conf_thresh <float>\n"
 64        "      Threshold for two images are from the same panorama confidence.\n"
 65        "      The default is 1.0.\n"
 66        "  --ba (reproj|ray)\n"
 67        "      Bundle adjustment cost function. The default is ray.\n"
 68        "  --ba_refine_mask (mask)\n"
 69        "      Set refinement mask for bundle adjustment. It looks like 'x_xxx',\n"
 70        "      where 'x' means refine respective parameter and '_' means don't\n"
 71        "      refine one, and has the following format:\n"
 72        "      <fx><skew><ppx><aspect><ppy>. The default mask is 'xxxxx'. If bundle\n"
 73        "      adjustment doesn't support estimation of selected parameter then\n"
 74        "      the respective flag is ignored.\n"
 75        "  --wave_correct (no|horiz|vert)\n"
 76        "      Perform wave effect correction. The default is 'horiz'.\n"
 77        "  --save_graph <file_name>\n"
 78        "      Save matches graph represented in DOT language to <file_name> file.\n"
 79        "      Labels description: Nm is number of matches, Ni is number of inliers,\n"
 80        "      C is confidence.\n"
 81        "\nCompositing Flags:\n"
 82        "  --warp (plane|cylindrical|spherical|fisheye|stereographic|compressedPlaneA2B1|compressedPlaneA1.5B1|compressedPlanePortraitA2B1|compressedPlanePortraitA1.5B1|paniniA2B1|paniniA1.5B1|paniniPortraitA2B1|paniniPortraitA1.5B1|mercator|transverseMercator)\n"
 83        "      Warp surface type. The default is 'spherical'.\n"
 84        "  --seam_megapix <float>\n"
 85        "      Resolution for seam estimation step. The default is 0.1 Mpx.\n"
 86        "  --seam (no|voronoi|gc_color|gc_colorgrad)\n"
 87        "      Seam estimation method. The default is 'gc_color'.\n"
 88        "  --compose_megapix <float>\n"
 89        "      Resolution for compositing step. Use -1 for original resolution.\n"
 90        "      The default is -1.\n"
 91        "  --expos_comp (no|gain|gain_blocks)\n"
 92        "      Exposure compensation method. The default is 'gain_blocks'.\n"
 93        "  --blend (no|feather|multiband)\n"
 94        "      Blending method. The default is 'multiband'.\n"
 95        "  --blend_strength <float>\n"
 96        "      Blending strength from [0,100] range. The default is 5.\n"
 97        "  --output <result_img>\n"
 98        "      The default is 'result.jpg'.\n";
 99}
100
101
102// Default command line args
103vector<string> img_names;
104bool preview = false;
105bool try_gpu = false;
106double work_megapix = 0.6;
107double seam_megapix = 0.1;
108double compose_megapix = -1;
109float conf_thresh = 1.f;
110string features = "surf";
111string ba_cost_func = "ray";
112string ba_refine_mask = "xxxxx";
113bool do_wave_correct = true;
114WaveCorrectKind wave_correct = detail::WAVE_CORRECT_HORIZ;
115bool save_graph = false;
116std::string save_graph_to;
117string warp_type = "spherical";
118int expos_comp_type = ExposureCompensator::GAIN_BLOCKS;
119float match_conf = 0.3f;
120string seam_find_type = "gc_color";
121int blend_type = Blender::MULTI_BAND;
122float blend_strength = 5;
123string result_name = "result.jpg";
124
125int parseCmdArgs(int argc, char** argv)
126{
127    if (argc == 1)
128    {
129        printUsage();
130        return -1;
131    }
132    for (int i = 1; i < argc; ++i)
133    {
134        if (string(argv[i]) == "--help" || string(argv[i]) == "/?")
135        {
136            printUsage();
137            return -1;
138        }
139        else if (string(argv[i]) == "--preview")
140        {
141            preview = true;
142        }
143        else if (string(argv[i]) == "--try_gpu")
144        {
145            if (string(argv[i + 1]) == "no")
146                try_gpu = false;
147            else if (string(argv[i + 1]) == "yes")
148                try_gpu = true;
149            else
150            {
151                cout << "Bad --try_gpu flag value\n";
152                return -1;
153            }
154            i++;
155        }
156        else if (string(argv[i]) == "--work_megapix")
157        {
158            work_megapix = atof(argv[i + 1]);
159            i++;
160        }
161        else if (string(argv[i]) == "--seam_megapix")
162        {
163            seam_megapix = atof(argv[i + 1]);
164            i++;
165        }
166        else if (string(argv[i]) == "--compose_megapix")
167        {
168            compose_megapix = atof(argv[i + 1]);
169            i++;
170        }
171        else if (string(argv[i]) == "--result")
172        {
173            result_name = argv[i + 1];
174            i++;
175        }
176        else if (string(argv[i]) == "--features")
177        {
178            features = argv[i + 1];
179            if (features == "orb")
180                match_conf = 0.3f;
181            i++;
182        }
183        else if (string(argv[i]) == "--match_conf")
184        {
185            match_conf = static_cast<float>(atof(argv[i + 1]));
186            i++;
187        }
188        else if (string(argv[i]) == "--conf_thresh")
189        {
190            conf_thresh = static_cast<float>(atof(argv[i + 1]));
191            i++;
192        }
193        else if (string(argv[i]) == "--ba")
194        {
195            ba_cost_func = argv[i + 1];
196            i++;
197        }
198        else if (string(argv[i]) == "--ba_refine_mask")
199        {
200            ba_refine_mask = argv[i + 1];
201            if (ba_refine_mask.size() != 5)
202            {
203                cout << "Incorrect refinement mask length.\n";
204                return -1;
205            }
206            i++;
207        }
208        else if (string(argv[i]) == "--wave_correct")
209        {
210            if (string(argv[i + 1]) == "no")
211                do_wave_correct = false;
212            else if (string(argv[i + 1]) == "horiz")
213            {
214                do_wave_correct = true;
215                wave_correct = detail::WAVE_CORRECT_HORIZ;
216            }
217            else if (string(argv[i + 1]) == "vert")
218            {
219                do_wave_correct = true;
220                wave_correct = detail::WAVE_CORRECT_VERT;
221            }
222            else
223            {
224                cout << "Bad --wave_correct flag value\n";
225                return -1;
226            }
227            i++;
228        }
229        else if (string(argv[i]) == "--save_graph")
230        {
231            save_graph = true;
232            save_graph_to = argv[i + 1];
233            i++;
234        }
235        else if (string(argv[i]) == "--warp")
236        {
237            warp_type = string(argv[i + 1]);
238            i++;
239        }
240        else if (string(argv[i]) == "--expos_comp")
241        {
242            if (string(argv[i + 1]) == "no")
243                expos_comp_type = ExposureCompensator::NO;
244            else if (string(argv[i + 1]) == "gain")
245                expos_comp_type = ExposureCompensator::GAIN;
246            else if (string(argv[i + 1]) == "gain_blocks")
247                expos_comp_type = ExposureCompensator::GAIN_BLOCKS;
248            else
249            {
250                cout << "Bad exposure compensation method\n";
251                return -1;
252            }
253            i++;
254        }
255        else if (string(argv[i]) == "--seam")
256        {
257            if (string(argv[i + 1]) == "no" ||
258                string(argv[i + 1]) == "voronoi" ||
259                string(argv[i + 1]) == "gc_color" ||
260                string(argv[i + 1]) == "gc_colorgrad")
261                seam_find_type = argv[i + 1];
262            else
263            {
264                cout << "Bad seam finding method\n";
265                return -1;
266            }
267            i++;
268        }
269        else if (string(argv[i]) == "--blend")
270        {
271            if (string(argv[i + 1]) == "no")
272                blend_type = Blender::NO;
273            else if (string(argv[i + 1]) == "feather")
274                blend_type = Blender::FEATHER;
275            else if (string(argv[i + 1]) == "multiband")
276                blend_type = Blender::MULTI_BAND;
277            else
278            {
279                cout << "Bad blending method\n";
280                return -1;
281            }
282            i++;
283        }
284        else if (string(argv[i]) == "--blend_strength")
285        {
286            blend_strength = static_cast<float>(atof(argv[i + 1]));
287            i++;
288        }
289        else if (string(argv[i]) == "--output")
290        {
291            result_name = argv[i + 1];
292            i++;
293        }
294        else
295            img_names.push_back(argv[i]);
296    }
297    if (preview)
298    {
299        compose_megapix = 0.6;
300    }
301    return 0;
302}
303
304
305int main(int argc, char* argv[])
306{
307    int64 app_start_time = getTickCount();
308    cv::setBreakOnError(true);
309
310    int retval = parseCmdArgs(argc, argv);
311    if (retval)
312        return retval;
313
314    // Check if have enough images
315    int num_images = static_cast<int>(img_names.size());
316    if (num_images < 2)
317    {
318        LOGLN("Need more images");
319        return -1;
320    }
321
322    double work_scale = 1, seam_scale = 1, compose_scale = 1;
323    bool is_work_scale_set = false, is_seam_scale_set = false, is_compose_scale_set = false;
324
325    LOGLN("Finding features...");
326    int64 t = getTickCount();
327
328    Ptr<FeaturesFinder> finder;
329    if (features == "surf")
330    {
331#ifdef HAVE_OPENCV_GPU
332        if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
333            finder = new SurfFeaturesFinderGpu();
334        else
335#endif
336            finder = new SurfFeaturesFinder();
337    }
338    else if (features == "orb")
339    {
340        finder = new OrbFeaturesFinder();
341    }
342    else
343    {
344        cout << "Unknown 2D features type: '" << features << "'.\n";
345        return -1;
346    }
347
348    Mat full_img, img;
349    vector<ImageFeatures> features(num_images);
350    vector<Mat> images(num_images);
351    vector<Size> full_img_sizes(num_images);
352    double seam_work_aspect = 1;
353
354    for (int i = 0; i < num_images; ++i)
355    {
356        full_img = imread(img_names[i]);
357        full_img_sizes[i] = full_img.size();
358
359        if (full_img.empty())
360        {
361            LOGLN("Can't open image " << img_names[i]);
362            return -1;
363        }
364        if (work_megapix < 0)
365        {
366            img = full_img;
367            work_scale = 1;
368            is_work_scale_set = true;
369        }
370        else
371        {
372            if (!is_work_scale_set)
373            {
374                work_scale = min(1.0, sqrt(work_megapix * 1e6 / full_img.size().area()));
375                is_work_scale_set = true;
376            }
377            resize(full_img, img, Size(), work_scale, work_scale);
378        }
379        if (!is_seam_scale_set)
380        {
381            seam_scale = min(1.0, sqrt(seam_megapix * 1e6 / full_img.size().area()));
382            seam_work_aspect = seam_scale / work_scale;
383            is_seam_scale_set = true;
384        }
385
386        (*finder)(img, features[i]);
387        features[i].img_idx = i;
388        LOGLN("Features in image #" << i+1 << ": " << features[i].keypoints.size());
389
390        resize(full_img, img, Size(), seam_scale, seam_scale);
391        images[i] = img.clone();
392    }
393
394    finder->collectGarbage();
395    full_img.release();
396    img.release();
397
398    LOGLN("Finding features, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
399
400    LOG("Pairwise matching");
401    t = getTickCount();
402    vector<MatchesInfo> pairwise_matches;
403    BestOf2NearestMatcher matcher(try_gpu, match_conf);
404    matcher(features, pairwise_matches);
405    matcher.collectGarbage();
406    LOGLN("Pairwise matching, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
407
408    // Check if we should save matches graph
409    if (save_graph)
410    {
411        LOGLN("Saving matches graph...");
412        ofstream f(save_graph_to.c_str());
413        f << matchesGraphAsString(img_names, pairwise_matches, conf_thresh);
414//        for (int i=0; i!=pairwise_matches.size(); ++i)
415//            f << pairwise_matches[i].src_img_idx << " " << pairwise_matches[i].dst_img_idx << "\n";
416    }
417
418    // Leave only images we are sure are from the same panorama
419    vector<int> indices = leaveBiggestComponent(features, pairwise_matches, conf_thresh);
420    vector<Mat> img_subset;
421    vector<string> img_names_subset;
422    vector<Size> full_img_sizes_subset;
423    for (size_t i = 0; i < indices.size(); ++i)
424    {
425        img_names_subset.push_back(img_names[indices[i]]);
426        img_subset.push_back(images[indices[i]]);
427        full_img_sizes_subset.push_back(full_img_sizes[indices[i]]);
428    }
429
430    images = img_subset;
431    img_names = img_names_subset;
432    full_img_sizes = full_img_sizes_subset;
433
434    // Check if we still have enough images
435    num_images = static_cast<int>(img_names.size());
436    if (num_images < 2)
437    {
438        LOGLN("Need more images");
439        return -1;
440    }
441
442    LOG("Homography-based init\n");
443    t = getTickCount();
444    HomographyBasedEstimator estimator;
445    vector<CameraParams> cameras;
446    estimator(features, pairwise_matches, cameras);
447    LOGLN("Homography-based init, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
448
449    for (size_t i = 0; i < cameras.size(); ++i)
450    {
451        Mat R;
452        cameras[i].R.convertTo(R, CV_32F);
453        cameras[i].R = R;
454        LOGLN("Initial intrinsics #" << indices[i]+1 << ":\n" << cameras[i].K());
455    }
456
457    LOG("Bundle Adjustment\n");
458    t = getTickCount();
459    Ptr<detail::BundleAdjusterBase> adjuster;
460    if (ba_cost_func == "reproj") adjuster = new detail::BundleAdjusterReproj();
461    else if (ba_cost_func == "ray") adjuster = new detail::BundleAdjusterRay();
462    else 
463    { 
464        cout << "Unknown bundle adjustment cost function: '" << ba_cost_func << "'.\n"; 
465        return -1; 
466    }
467    adjuster->setConfThresh(conf_thresh);
468    Mat_<uchar> refine_mask = Mat::zeros(3, 3, CV_8U);
469    if (ba_refine_mask[0] == 'x') refine_mask(0,0) = 1;
470    if (ba_refine_mask[1] == 'x') refine_mask(0,1) = 1;
471    if (ba_refine_mask[2] == 'x') refine_mask(0,2) = 1;
472    if (ba_refine_mask[3] == 'x') refine_mask(1,1) = 1;
473    if (ba_refine_mask[4] == 'x') refine_mask(1,2) = 1;
474    adjuster->setRefinementMask(refine_mask);
475    (*adjuster)(features, pairwise_matches, cameras);
476    LOGLN("Bundle Adjustment, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
477
478    // Find median focal length
479
480    vector<double> focals;
481    for (size_t i = 0; i < cameras.size(); ++i)
482    {
483        LOGLN("Camera #" << indices[i]+1 << ":\n" << cameras[i].K());
484        focals.push_back(cameras[i].focal);
485    }
486
487    sort(focals.begin(), focals.end());
488    float warped_image_scale;
489    if (focals.size() % 2 == 1)
490        warped_image_scale = static_cast<float>(focals[focals.size() / 2]);
491    else
492        warped_image_scale = static_cast<float>(focals[focals.size() / 2 - 1] + focals[focals.size() / 2]) * 0.5f;
493
494    if (do_wave_correct)
495    {
496        vector<Mat> rmats;
497        for (size_t i = 0; i < cameras.size(); ++i)
498            rmats.push_back(cameras[i].R);
499        waveCorrect(rmats, wave_correct);
500        for (size_t i = 0; i < cameras.size(); ++i)
501            cameras[i].R = rmats[i];
502    }
503
504    LOGLN("Warping images (auxiliary)... ");
505    t = getTickCount();
506
507    vector<Point> corners(num_images);
508    vector<Mat> masks_warped(num_images);
509    vector<Mat> images_warped(num_images);
510    vector<Size> sizes(num_images);
511    vector<Mat> masks(num_images);
512
513    // Preapre images masks
514    for (int i = 0; i < num_images; ++i)
515    {
516        masks[i].create(images[i].size(), CV_8U);
517        masks[i].setTo(Scalar::all(255));
518    }
519
520    // Warp images and their masks
521
522    Ptr<WarperCreator> warper_creator;
523#ifdef HAVE_OPENCV_GPU
524    if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
525    {
526        if (warp_type == "plane") warper_creator = new cv::PlaneWarperGpu();
527        else if (warp_type == "cylindrical") warper_creator = new cv::CylindricalWarperGpu();
528        else if (warp_type == "spherical") warper_creator = new cv::SphericalWarperGpu();
529    }
530    else
531#endif
532    {
533        if (warp_type == "plane") warper_creator = new cv::PlaneWarper();
534        else if (warp_type == "cylindrical") warper_creator = new cv::CylindricalWarper();
535        else if (warp_type == "spherical") warper_creator = new cv::SphericalWarper();
536		else if (warp_type == "fisheye") warper_creator = new cv::FisheyeWarper();
537		else if (warp_type == "stereographic") warper_creator = new cv::StereographicWarper();
538		else if (warp_type == "compressedPlaneA2B1") warper_creator = new cv::CompressedRectilinearWarper(2, 1);
539		else if (warp_type == "compressedPlaneA1.5B1") warper_creator = new cv::CompressedRectilinearWarper(1.5, 1);
540		else if (warp_type == "compressedPlanePortraitA2B1") warper_creator = new cv::CompressedRectilinearPortraitWarper(2, 1);
541		else if (warp_type == "compressedPlanePortraitA1.5B1") warper_creator = new cv::CompressedRectilinearPortraitWarper(1.5, 1);
542		else if (warp_type == "paniniA2B1") warper_creator = new cv::PaniniWarper(2, 1);
543		else if (warp_type == "paniniA1.5B1") warper_creator = new cv::PaniniWarper(1.5, 1);
544		else if (warp_type == "paniniPortraitA2B1") warper_creator = new cv::PaniniPortraitWarper(2, 1);
545		else if (warp_type == "paniniPortraitA1.5B1") warper_creator = new cv::PaniniPortraitWarper(1.5, 1);
546		else if (warp_type == "mercator") warper_creator = new cv::MercatorWarper();
547		else if (warp_type == "transverseMercator") warper_creator = new cv::TransverseMercatorWarper();
548    }
549
550    if (warper_creator.empty())
551    {
552        cout << "Can't create the following warper '" << warp_type << "'\n";
553        return 1;
554    }
555    
556    Ptr<RotationWarper> warper = warper_creator->create(static_cast<float>(warped_image_scale * seam_work_aspect));
557
558    for (int i = 0; i < num_images; ++i)
559    {
560        Mat_<float> K;
561        cameras[i].K().convertTo(K, CV_32F);
562        float swa = (float)seam_work_aspect;
563        K(0,0) *= swa; K(0,2) *= swa;
564        K(1,1) *= swa; K(1,2) *= swa;
565
566        corners[i] = warper->warp(images[i], K, cameras[i].R, INTER_LINEAR, BORDER_REFLECT, images_warped[i]);
567        sizes[i] = images_warped[i].size();
568
569        warper->warp(masks[i], K, cameras[i].R, INTER_NEAREST, BORDER_CONSTANT, masks_warped[i]);
570    }
571
572    vector<Mat> images_warped_f(num_images);
573    for (int i = 0; i < num_images; ++i)
574        images_warped[i].convertTo(images_warped_f[i], CV_32F);
575
576    LOGLN("Warping images, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
577
578    Ptr<ExposureCompensator> compensator = ExposureCompensator::createDefault(expos_comp_type);
579    compensator->feed(corners, images_warped, masks_warped);
580
581    Ptr<SeamFinder> seam_finder;
582    if (seam_find_type == "no")
583        seam_finder = new detail::NoSeamFinder();
584    else if (seam_find_type == "voronoi")
585        seam_finder = new detail::VoronoiSeamFinder();
586    else if (seam_find_type == "gc_color")
587    {
588#ifdef HAVE_OPENCV_GPU
589        if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
590            seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR);
591        else
592#endif
593            seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR);
594    }
595    else if (seam_find_type == "gc_colorgrad")
596    {
597#ifdef HAVE_OPENCV_GPU
598        if (try_gpu && gpu::getCudaEnabledDeviceCount() > 0)
599            seam_finder = new detail::GraphCutSeamFinderGpu(GraphCutSeamFinderBase::COST_COLOR_GRAD);
600        else
601#endif
602            seam_finder = new detail::GraphCutSeamFinder(GraphCutSeamFinderBase::COST_COLOR_GRAD);
603    }
604    if (seam_finder.empty())
605    {
606        cout << "Can't create the following seam finder '" << seam_find_type << "'\n";
607        return 1;
608    }
609
610    seam_finder->find(images_warped_f, corners, masks_warped);
611
612    // Release unused memory
613    images.clear();
614    images_warped.clear();
615    images_warped_f.clear();
616    masks.clear();
617
618    LOGLN("Compositing...");
619    t = getTickCount();
620
621    Mat img_warped, img_warped_s;
622    Mat dilated_mask, seam_mask, mask, mask_warped;
623    Ptr<Blender> blender;
624    double compose_seam_aspect = 1;
625    double compose_work_aspect = 1;
626
627    for (int img_idx = 0; img_idx < num_images; ++img_idx)
628    {
629        LOGLN("Compositing image #" << indices[img_idx]+1);
630
631        // Read image and resize it if necessary
632        full_img = imread(img_names[img_idx]);
633        if (!is_compose_scale_set)
634        {
635            if (compose_megapix > 0)
636                compose_scale = min(1.0, sqrt(compose_megapix * 1e6 / full_img.size().area()));
637            is_compose_scale_set = true;
638
639            // Compute relative scales
640            compose_seam_aspect = compose_scale / seam_scale;
641            compose_work_aspect = compose_scale / work_scale;
642
643            // Update warped image scale
644            warped_image_scale *= static_cast<float>(compose_work_aspect);
645            warper = warper_creator->create(warped_image_scale);
646
647            // Update corners and sizes
648            for (int i = 0; i < num_images; ++i)
649            {
650                // Update intrinsics
651                cameras[i].focal *= compose_work_aspect;
652                cameras[i].ppx *= compose_work_aspect;
653                cameras[i].ppy *= compose_work_aspect;
654
655                // Update corner and size
656                Size sz = full_img_sizes[i];
657                if (std::abs(compose_scale - 1) > 1e-1)
658                {
659                    sz.width = cvRound(full_img_sizes[i].width * compose_scale);
660                    sz.height = cvRound(full_img_sizes[i].height * compose_scale);
661                }
662
663                Mat K;
664                cameras[i].K().convertTo(K, CV_32F);
665                Rect roi = warper->warpRoi(sz, K, cameras[i].R);
666                corners[i] = roi.tl();
667                sizes[i] = roi.size();
668            }
669        }
670        if (abs(compose_scale - 1) > 1e-1)
671            resize(full_img, img, Size(), compose_scale, compose_scale);
672        else
673            img = full_img;
674        full_img.release();
675        Size img_size = img.size();
676
677        Mat K;
678        cameras[img_idx].K().convertTo(K, CV_32F);
679
680        // Warp the current image
681        warper->warp(img, K, cameras[img_idx].R, INTER_LINEAR, BORDER_REFLECT, img_warped);
682
683        // Warp the current image mask
684        mask.create(img_size, CV_8U);
685        mask.setTo(Scalar::all(255));
686        warper->warp(mask, K, cameras[img_idx].R, INTER_NEAREST, BORDER_CONSTANT, mask_warped);
687
688        // Compensate exposure
689        compensator->apply(img_idx, corners[img_idx], img_warped, mask_warped);
690
691        img_warped.convertTo(img_warped_s, CV_16S);
692        img_warped.release();
693        img.release();
694        mask.release();
695
696        dilate(masks_warped[img_idx], dilated_mask, Mat());
697        resize(dilated_mask, seam_mask, mask_warped.size());
698        mask_warped = seam_mask & mask_warped;
699
700        if (blender.empty())
701        {
702            blender = Blender::createDefault(blend_type, try_gpu);
703            Size dst_sz = resultRoi(corners, sizes).size();
704            float blend_width = sqrt(static_cast<float>(dst_sz.area())) * blend_strength / 100.f;
705            if (blend_width < 1.f)
706                blender = Blender::createDefault(Blender::NO, try_gpu);
707            else if (blend_type == Blender::MULTI_BAND)
708            {
709                MultiBandBlender* mb = dynamic_cast<MultiBandBlender*>(static_cast<Blender*>(blender));
710                mb->setNumBands(static_cast<int>(ceil(log(blend_width)/log(2.)) - 1.));
711                LOGLN("Multi-band blender, number of bands: " << mb->numBands());
712            }
713            else if (blend_type == Blender::FEATHER)
714            {
715                FeatherBlender* fb = dynamic_cast<FeatherBlender*>(static_cast<Blender*>(blender));
716                fb->setSharpness(1.f/blend_width);
717                LOGLN("Feather blender, sharpness: " << fb->sharpness());
718            }
719            blender->prepare(corners, sizes);
720        }
721
722        // Blend the current image
723        blender->feed(img_warped_s, mask_warped, corners[img_idx]);
724    }
725
726    Mat result, result_mask;
727    blender->blend(result, result_mask);
728
729    LOGLN("Compositing, time: " << ((getTickCount() - t) / getTickFrequency()) << " sec");
730
731    imwrite(result_name, result);
732
733    LOGLN("Finished, total time: " << ((getTickCount() - app_start_time) / getTickFrequency()) << " sec");
734    return 0;
735}
736
737