align_image_stack.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00002 
00027 #include <hugin_config.h>
00028 #include <fstream>
00029 #include <sstream>
00030 #include <iostream>
00031 
00032 #include <vigra/error.hxx>
00033 #include <vigra/cornerdetection.hxx>
00034 #include <vigra/localminmax.hxx>
00035 #include <hugin_utils/utils.h>
00036 
00037 #include <vigra_ext/Pyramid.h>
00038 #include <vigra_ext/Correlation.h>
00039 #include <vigra_ext/InterestPoints.h>
00040 #include <vigra_ext/impexalpha.hxx>
00041 
00042 #include <panodata/Panorama.h>
00043 #include <panotools/PanoToolsOptimizerWrapper.h>
00044 #include <panodata/StandardImageVariableGroups.h>
00045 #include <algorithms/optimizer/PTOptimizer.h>
00046 #include <nona/Stitcher.h>
00047 #include <algorithms/basic/CalculateOptimalROI.h>
00048 #include <lensdb/LensDB.h>
00049 
00050 #include <getopt.h>
00051 
00052 #include <hugin_utils/openmp_lock.h>
00053 #include <tiff.h>
00054 
00055 #ifdef __APPLE__
00056 #include <hugin_config.h>
00057 #include <mach-o/dyld.h>    /* _NSGetExecutablePath */
00058 #include <limits.h>         /* PATH_MAX */
00059 #include <libgen.h>         /* dirname */
00060 #endif
00061 
00062 int g_verbose = 0;
00063 
00064 static void usage(const char* name)
00065 {
00066     std::cerr << name << ": align overlapping images for HDR creation" << std::endl
00067          << "align_image_stack version " << hugin_utils::GetHuginVersion() << std::endl
00068          << std::endl
00069          << "Usage: " << name  << " [options] input files" << std::endl
00070          << "Valid options are:" << std::endl
00071          << " Modes of operation:" << std::endl
00072          << "  -p file   Output .pto file (useful for debugging, or further refinement)" << std::endl
00073          << "  -a prefix align images, output as prefix_xxxx.tif" << std::endl
00074          << "  -o output merge images to HDR, generate output.hdr" << std::endl
00075          << " Modifiers" << std::endl
00076          << "  -v        Verbose, print progress messages.  Repeat for higher verbosity" << std::endl
00077          << "  -e        Assume input images are full frame fish eye (default: rectilinear)" << std::endl
00078          << "  -t num    Remove all control points with an error higher than num pixels" << std::endl
00079          << "             (default: 3)" << std::endl
00080          << "  --corr=num  Correlation threshold for identifing control points" << std::endl
00081          << "               (default: 0.9)" << std::endl
00082          << "  -f HFOV   approximate horizontal field of view of input images," << std::endl
00083          << "             use if EXIF info not complete" << std::endl
00084          << "  -m        Optimize field of view for all images, except for first." << std::endl
00085          << "             Useful for aligning focus stacks with slightly" << std::endl
00086          << "             different magnification." << std::endl
00087          << "  -d        Optimize radial distortion for all images, except for first." << std::endl
00088          << "  -i        Optimize image center shift for all images, except for first." << std::endl
00089          << "  -x        Optimize X coordinate of the camera position." << std::endl
00090          << "  -y        Optimize Y coordinate of the camera position." << std::endl
00091          << "  -z        Optimize Z coordinate of the camera position." << std::endl
00092          << "             Useful for aligning more distorted images." << std::endl
00093          << "  -S        Assume stereo images - allow horizontal shift of control points." << std::endl
00094          << "  -A        Align stereo window - assumes -S." << std::endl
00095          << "  -P        Align stereo window with pop-out effect - assumes -S." << std::endl
00096          << "  -C        Auto crop the image to the area covered by all images." << std::endl
00097          << "  -c num    number of control points (per grid) to create" << std::endl
00098          << "             between adjacent images (default: 8)" << std::endl
00099          << "  -l        Assume linear input files" << std::endl
00100          << "  -s scale  Scale down image by 2^scale (default: 1 [2x downsampling])" << std::endl
00101          << "  -g gsize  Break image into a rectangular grid (gsize x gsize) and attempt" << std::endl
00102          << "             to find num control points in each section" << std::endl
00103          << "             (default: 5 [5x5 grid] )" << std::endl
00104          << "  --distortion Try to load distortion information from lens database" << std::endl
00105          << "  --use-given-order  Use the image order as given from command line" << std::endl
00106          << "                     (By default images will be sorted by exposure values.)" << std::endl
00107          << "  --gpu     Use GPU for remapping" << std::endl
00108          << "  -h        Display help (this text)" << std::endl
00109          << std::endl;
00110 }
00111 
00112 typedef std::multimap<double, vigra::Diff2D> MapPoints;
00113 static hugin_omp::Lock lock;
00114 
00115 namespace detail
00116 {
00117     template <class ImageType>
00118     vigra_ext::CorrelationResult FineTunePoint(const ImageType& leftImg, const vigra::Diff2D templPos, const int templSize,
00119         const ImageType& rightImg, const vigra::Diff2D searchPos, const int searchWidth, vigra::VigraTrueType)
00120     {
00121         return vigra_ext::PointFineTune(leftImg, leftImg.accessor(),
00122             templPos, templSize,
00123             rightImg, rightImg.accessor(),
00124             searchPos, searchWidth);
00125     };
00126 
00127     template <class ImageType>
00128     vigra_ext::CorrelationResult FineTunePoint(const ImageType& leftImg, const vigra::Diff2D templPos, const int templSize,
00129         const ImageType& rightImg, const vigra::Diff2D searchPos, const int searchWidth, vigra::VigraFalseType)
00130     {
00131         return vigra_ext::PointFineTune(leftImg,
00132             vigra::RGBToGrayAccessor<typename ImageType::value_type>(),
00133             templPos, templSize,
00134             rightImg, vigra::RGBToGrayAccessor<typename ImageType::value_type>(),
00135             searchPos, searchWidth);
00136     };
00137 }
00138 template <class ImageType>
00139 void FineTuneInterestPoints(HuginBase::Panorama& pano,
00140                             int img1, const ImageType& leftImg, const ImageType& leftImgOrig,
00141                             int img2, const ImageType& rightImg, const ImageType& rightImgOrig,
00142                             const MapPoints& points, unsigned nPoints, int pyrLevel, int templWidth, int sWidth, double scaleFactor, double corrThresh, bool stereo)
00143 {
00144     typedef typename ImageType::value_type ImageValueType;
00145     typedef typename vigra::NumericTraits<ImageValueType>::isScalar is_scalar;
00146 
00147     unsigned nGood = 0;
00148     // loop over all points, starting with the highest corner score
00149     for (MapPoints::const_reverse_iterator it = points.rbegin(); it != points.rend(); ++it)
00150     {
00151         if (nGood >= nPoints)
00152         {
00153             // we have enough points, stop
00154             break;
00155         }
00156 
00157         vigra_ext::CorrelationResult res = detail::FineTunePoint(leftImg, it->second, templWidth,
00158             rightImg, it->second, sWidth, is_scalar());
00159         if (g_verbose > 2)
00160         {
00161             std::ostringstream buf;
00162             buf << "I :" << (*it).second.x* scaleFactor << "," << (*it).second.y* scaleFactor << " -> "
00163                 << res.maxpos.x* scaleFactor << "," << res.maxpos.y* scaleFactor << ":  corr coeff: " << res.maxi
00164                 << " curv:" << res.curv.x << " " << res.curv.y << std::endl;
00165             std::cout << buf.str();
00166         }
00167         if (res.maxi < corrThresh)
00168         {
00169             DEBUG_DEBUG("low correlation: " << res.maxi << " curv: " << res.curv);
00170             continue;
00171         }
00172 
00173         if (pyrLevel > 0)
00174         {
00175             res = detail::FineTunePoint(leftImgOrig, vigra::Diff2D((*it).second.x * scaleFactor, (*it).second.y * scaleFactor),
00176                 templWidth, rightImgOrig, vigra::Diff2D(res.maxpos.x * scaleFactor, res.maxpos.y * scaleFactor),
00177                 scaleFactor, is_scalar());
00178 
00179             if (g_verbose > 2)
00180             {
00181                 std::ostringstream buf;
00182                 buf << "II>" << (*it).second.x* scaleFactor << "," << (*it).second.y* scaleFactor << " -> "
00183                     << res.maxpos.x << "," << res.maxpos.y << ":  corr coeff: " << res.maxi
00184                     << " curv:" << res.curv.x << " " << res.curv.y << std::endl;
00185                 std::cout << buf.str();
00186             }
00187             if (res.maxi < corrThresh)
00188             {
00189                 DEBUG_DEBUG("low correlation in pass 2: " << res.maxi << " curv: " << res.curv);
00190                 continue;
00191             }
00192         }
00193 
00194         nGood++;
00195         // add control point
00196         HuginBase::ControlPoint p(img1, (*it).second.x * scaleFactor,
00197                        (*it).second.y * scaleFactor,
00198                        img2, res.maxpos.x,
00199                        res.maxpos.y,
00200                        stereo ? HuginBase::ControlPoint::Y : HuginBase::ControlPoint::X_Y);
00201         {
00202             hugin_omp::ScopedLock sl(lock);
00203             pano.addCtrlPoint(p);
00204         };
00205     }
00206     if (g_verbose > 0)
00207     {
00208         std::ostringstream buf;
00209         buf << "Number of good matches: " << nGood << ", bad matches: " << points.size() - nGood << std::endl;
00210         std::cout << buf.str();
00211     }
00212 };
00213 
00214 namespace detail
00215 {
00216     template <class ImageType>
00217     void FindInterestPointsPartial(const ImageType& image, const vigra::Rect2D& rect, double scale,
00218         unsigned nPoints, std::multimap<double, vigra::Diff2D> &points, vigra::VigraTrueType)
00219     {
00220         vigra_ext::findInterestPointsPartial(vigra::srcImageRange(image), rect, scale, nPoints, points);
00221     };
00222 
00223     template <class ImageType>
00224     void FindInterestPointsPartial(const ImageType& image, const vigra::Rect2D& rect, double scale,
00225         unsigned nPoints, std::multimap<double, vigra::Diff2D> &points, vigra::VigraFalseType)
00226     {
00227         typedef typename ImageType::value_type ImageValueType;
00228         vigra_ext::findInterestPointsPartial(vigra::srcImageRange(image, vigra::RGBToGrayAccessor<ImageValueType>()), rect, scale, nPoints, points);
00229     };
00230 };
00231 
00232 template <class ImageType>
00233 void createCtrlPoints(HuginBase::Panorama& pano, int img1, const ImageType& leftImg, const ImageType& leftImgOrig, int img2, const ImageType& rightImg, const ImageType& rightImgOrig, int pyrLevel, double scale, unsigned nPoints, unsigned grid, double corrThresh = 0.9, bool stereo = false)
00234 {
00235     typedef typename ImageType::value_type ImageValueType;
00236     typedef typename vigra::NumericTraits<ImageValueType>::isScalar is_scalar;
00237 
00239     // find interesting corners using harris corner detector
00240     if (stereo)
00241     {
00242         // add one vertical control point to keep the images aligned vertically
00243         HuginBase::ControlPoint p(img1, 0, 0, img2, 0, 0, HuginBase::ControlPoint::X);
00244         pano.addCtrlPoint(p);
00245     }
00246 
00247     if (g_verbose > 0)
00248     {
00249         std::cout << "Trying to find " << nPoints << " corners... " << std::endl;
00250     }
00251 
00252     vigra::Size2D size(leftImg.width(), leftImg.height());
00253     std::vector<vigra::Rect2D> rects;
00254     for (unsigned party = 0; party < grid; party++)
00255     {
00256         for (unsigned partx = 0; partx < grid; partx++)
00257         {
00258             // run corner detector only in current sub-region (saves a lot of memory for big images)
00259             vigra::Rect2D rect(partx*size.x / grid, party*size.y / grid,
00260                                (partx + 1)*size.x / grid, (party + 1)*size.y / grid);
00261             rect &= vigra::Rect2D(size);
00262             if (rect.width()>0 && rect.height()>0)
00263             {
00264                 rects.push_back(rect);
00265             };
00266         };
00267     };
00268 
00269     const double scaleFactor = 1 << pyrLevel;
00270     const long templWidth = 20;
00271     const long sWidth = 100;
00272 
00273     #pragma omp parallel for schedule(dynamic)
00274     for (int i = 0; i < rects.size(); ++i)
00275     {
00276         MapPoints points;
00277         vigra::Rect2D rect(rects[i]);
00278         detail::FindInterestPointsPartial(leftImg, rect, scale, 5 * nPoints, points, is_scalar());
00279         FineTuneInterestPoints(pano, img1, leftImg, leftImgOrig, img2, rightImg, rightImgOrig, points, nPoints, pyrLevel, templWidth, sWidth, scaleFactor, corrThresh, stereo);
00280     };
00281 
00282     if (stereo)
00283     {
00284         // add some additional control points around image edges
00285         // this is useful for better results - images are more distorted around edges
00286         // and also for stereoscopic window adjustment - it must be alligned according to
00287         // the nearest object which crosses the edge and these control points helps to find it.
00288         MapPoints up;
00289         MapPoints down;
00290         MapPoints left;
00291         MapPoints right;
00292         int xstep = leftImg.size().x / (nPoints + 1);
00293         int ystep = leftImg.size().y / (nPoints + 1);
00294         for (int k = 6; k >= 0; --k)
00295         {
00296             for (int j = 0; j < 2; ++j)
00297             {
00298                 for (unsigned int i = 0; i < nPoints; ++i)
00299                 {
00300                     up.insert(std::make_pair(0, vigra::Diff2D(j * xstep / 2 + i * xstep, 1 + k * 10)));
00301                     down.insert(std::make_pair(0, vigra::Diff2D(j * xstep / 2 + i * xstep, leftImg.size().y - 2 - k * 10)));
00302                     left.insert(std::make_pair(0, vigra::Diff2D(1 + k * 10, j * ystep / 2 + i * ystep)));
00303                     right.insert(std::make_pair(0, vigra::Diff2D(leftImg.size().x - 2 - k * 10, j * ystep / 2 + i * ystep)));
00304                 };
00305             };
00306         };
00307         // execute parallel
00308         #pragma omp parallel sections
00309         {
00310             #pragma omp section
00311             {
00312                 FineTuneInterestPoints(pano, img1, leftImg, leftImgOrig, img2, rightImg, rightImgOrig, up, nPoints, pyrLevel, templWidth, sWidth, corrThresh, scaleFactor, stereo);
00313             }
00314             #pragma omp section
00315             {
00316                 FineTuneInterestPoints(pano, img1, leftImg, leftImgOrig, img2, rightImg, rightImgOrig, down, nPoints, pyrLevel, templWidth, sWidth, corrThresh, scaleFactor, stereo);
00317             }
00318             #pragma omp section
00319             {
00320                 FineTuneInterestPoints(pano, img1, leftImg, leftImgOrig, img2, rightImg, rightImgOrig, left, nPoints, pyrLevel, templWidth, sWidth, corrThresh, scaleFactor, stereo);
00321             }
00322             #pragma omp section
00323             {
00324                 FineTuneInterestPoints(pano, img1, leftImg, leftImgOrig, img2, rightImg, rightImgOrig, right, nPoints, pyrLevel, templWidth, sWidth, corrThresh, scaleFactor, stereo);
00325             }
00326         }
00327     }
00328 };
00329 
00330 void alignStereoWindow(HuginBase::Panorama& pano, bool pop_out)
00331 {
00332     HuginBase::CPVector cps = pano.getCtrlPoints();
00333     std::vector<HuginBase::PTools::Transform*> transTable(pano.getNrOfImages());
00334 
00335     std::vector<int> max_i(pano.getNrOfImages() - 1, -1); // index of a point with biggest shift
00336     std::vector<int> max_i_b(pano.getNrOfImages() - 1, -1); // the same as above, only border points considered
00337     std::vector<double> max_dif(pano.getNrOfImages() - 1, -1000000000); // value of the shift
00338     std::vector<double> max_dif_b(pano.getNrOfImages() - 1, -1000000000); // the same as above, only border points considered
00339 
00340     for (int i=0; i < pano.getNrOfImages(); i++)
00341     {
00342         transTable[i] = new HuginBase::PTools::Transform();
00343         transTable[i]->createInvTransform(pano.getImage(i), pano.getOptions());
00344     }
00345 
00346     double rbs = 0.1; // relative border size
00347 
00348     for (int i=0; i < (int)cps.size(); i++)
00349     {
00350         if (cps[i].mode == HuginBase::ControlPoint::X)
00351         {
00352             if (max_i[cps[i].image1Nr] < 0) // first control point for given pair
00353             {
00354                 max_i[cps[i].image1Nr] = i;    // use it as a fallback in case on other points exist
00355             }
00356             continue;
00357         }
00358 
00359         vigra::Size2D size1 = pano.getImage(cps[i].image1Nr).getSize();
00360         vigra::Size2D size2 = pano.getImage(cps[i].image2Nr).getSize();
00361 
00362         vigra::Rect2D rect1(size1);
00363         vigra::Rect2D rect2(size2);
00364 
00365         rect1.addBorder(-size1.width() * rbs, -size1.height() * rbs);
00366         rect2.addBorder(-size2.width() * rbs, -size2.height() * rbs);
00367 
00368 
00369         double xt1, yt1, xt2, yt2;
00370         if(!transTable[cps[i].image1Nr]->transformImgCoord(xt1, yt1, cps[i].x1, cps[i].y1))
00371         {
00372             continue;
00373         }
00374         if(!transTable[cps[i].image2Nr]->transformImgCoord(xt2, yt2, cps[i].x2, cps[i].y2))
00375         {
00376             continue;
00377         }
00378 
00379         double dif = xt2 - xt1;
00380         if (dif > max_dif[cps[i].image1Nr])
00381         {
00382             max_dif[cps[i].image1Nr] = dif;
00383             max_i[cps[i].image1Nr] = i;
00384         }
00385 
00386         if (!(rect1.contains(vigra::Point2D(cps[i].x1, cps[i].y1)) &&
00387                 rect2.contains(vigra::Point2D(cps[i].x2, cps[i].y2))))
00388         {
00389             // the same for border points only
00390             if (dif > max_dif_b[cps[i].image1Nr])
00391             {
00392                 max_dif_b[cps[i].image1Nr] = dif;
00393                 max_i_b[cps[i].image1Nr] = i;
00394             }
00395         }
00396     }
00397 
00398     for (int i=0; i < pano.getNrOfImages(); i++)
00399     {
00400         delete transTable[i];
00401     }
00402 
00403     for (int i=0; i < (int)max_i.size(); i++)
00404     {
00405         if (pop_out && (max_i_b[i] >= 0)) // check points at border
00406         {
00407             cps[max_i_b[i]].mode = HuginBase::ControlPoint::X_Y;
00408         }
00409         else if (max_i[i] >= 0) // no points at border - use any point
00410         {
00411             cps[max_i[i]].mode = HuginBase::ControlPoint::X_Y;
00412         }
00413         else
00414         {
00415             //no points at all - should not happen
00416         }
00417 
00418     }
00419 
00420     HuginBase::CPVector newCPs;
00421     for (int i=0; i < (int)cps.size(); i++)
00422     {
00423         if (cps[i].mode != HuginBase::ControlPoint::X)   // remove the vertical control lines, X_Y points replaces them
00424         {
00425             newCPs.push_back(cps[i]);
00426         }
00427     }
00428 
00429     pano.setCtrlPoints(newCPs);
00430 }
00431 
00432 void autoCrop(HuginBase::Panorama& pano)
00433 {
00434     AppBase::DummyProgressDisplay dummy;
00435     HuginBase::CalculateOptimalROI cropPano(pano, &dummy, true);
00436     cropPano.run();
00437 
00438     vigra::Rect2D roi=cropPano.getResultOptimalROI();
00439     //set the ROI - fail if the right/bottom is zero, meaning all zero
00440     if(!roi.isEmpty())
00441     {
00442         HuginBase::PanoramaOptions opt = pano.getOptions();
00443         opt.setROI(roi);
00444         pano.setOptions(opt);
00445         std::cout << "Set crop size to " << roi.left() << "," << roi.top() << "," << roi.right() << "," << roi.bottom() << std::endl;
00446     }
00447     else
00448     {
00449         std::cout << "Could not find best crop rectangle for image" << std::endl;
00450     };
00451 }
00452 
00453 struct Parameters
00454 {
00455     Parameters()
00456     {
00457         cpErrorThreshold = 3;
00458         corrThresh = 0.9;
00459         nPoints = 8;
00460         grid = 5;
00461         hfov = 0;
00462         pyrLevel = 1;
00463         linear = false;   // Assume linear input files if true
00464         optHFOV = false;
00465         optDistortion = false;
00466         optCenter = false;
00467         optX = false;
00468         optY = false;
00469         optZ = false;
00470         stereo = false;
00471         stereo_window = false;
00472         pop_out = false;
00473         crop = false;
00474         fisheye = false;
00475         gpu = false;
00476         loadDistortion = false;
00477         sortImagesByEv = true;
00478     }
00479 
00480     double cpErrorThreshold;
00481     double corrThresh;
00482     int nPoints;
00483     int grid;           // Partition images into grid x grid subregions, each with npoints
00484     double hfov;
00485     bool linear;
00486     bool optHFOV;
00487     bool optDistortion;
00488     bool optCenter;
00489     bool optX;
00490     bool optY;
00491     bool optZ;
00492     bool fisheye;
00493     bool stereo;
00494     bool stereo_window;
00495     bool pop_out;
00496     bool crop;
00497     bool gpu;
00498     bool loadDistortion;
00499     bool sortImagesByEv;
00500     int pyrLevel;
00501     std::string alignedPrefix;
00502     std::string ptoFile;
00503     std::string hdrFile;
00504     std::string basename;
00505 };
00506 
00507 // dummy panotools progress functions
00508 static int ptProgress(int command, char* argument)
00509 {
00510     return 1;
00511 }
00512 static int ptinfoDlg(int command, char* argument)
00513 {
00514     return 1;
00515 }
00516 
00517 struct SortImageVectorEV
00518 {
00519 public:
00520     explicit SortImageVectorEV(const HuginBase::Panorama* pano) : m_pano(pano) {};
00521     bool operator()(const unsigned int i, const unsigned int j)
00522     {
00523         return m_pano->getImage(i).getExposureValue() > m_pano->getImage(j).getExposureValue();
00524     };
00525 private:
00526     const HuginBase::Panorama* m_pano;
00527 };
00528 
00529 template <class PixelType>
00530 int main2(std::vector<std::string> files, Parameters param)
00531 {
00532     typedef vigra::BasicImage<PixelType> ImageType;
00533     try
00534     {
00535         HuginBase::Panorama pano;
00536         HuginBase::Lens l;
00537 
00538         // add the first image.to the panorama object
00539         HuginBase::SrcPanoImage srcImg;
00540         srcImg.setFilename(files[0]);
00541 
00542         if (param.fisheye)
00543         {
00544             srcImg.setProjection(HuginBase::SrcPanoImage::FULL_FRAME_FISHEYE);
00545         }
00546         srcImg.readEXIF();
00547         srcImg.applyEXIFValues();
00548         if (param.sortImagesByEv)
00549         {
00550             if (fabs(srcImg.getExposureValue()) < 1E-6)
00551             {
00552                 // no exposure values found in file, don't sort images
00553                 param.sortImagesByEv = false;
00554             }
00555         };
00556         // disable autorotate
00557         srcImg.setRoll(0);
00558         if (srcImg.getSize().x == 0 || srcImg.getSize().y == 0)
00559         {
00560             std::cerr << "Could not decode image: " << files[0] << "Unsupported image file format" << std::endl;
00561             return 1;
00562         }
00563 
00564         if(param.loadDistortion)
00565         {
00566             if(srcImg.readDistortionFromDB())
00567             {
00568                 std::cout << "\tRead distortion data from lens database." << std::endl;
00569             }
00570             else
00571             {
00572                 std::cout << "\tNo valid distortion data found in lens database." << std::endl;
00573             }
00574         }
00575 
00576         // use hfov specified by user.
00577         if (param.hfov > 0)
00578         {
00579             srcImg.setHFOV(param.hfov);
00580         }
00581         else if (srcImg.getCropFactor() == 0)
00582         {
00583             // could not read HFOV, assuming default: 50
00584             srcImg.setHFOV(50);
00585         }
00586 
00587         if (param.linear)
00588         {
00589             srcImg.setResponseType(HuginBase::SrcPanoImage::RESPONSE_LINEAR);
00590             if (g_verbose>0)
00591             {
00592                 std::cout << "Using linear response" << std::endl;
00593             }
00594         }
00595 
00596         pano.addImage(srcImg);
00597 
00598         // setup output to be exactly similar to input image
00599         HuginBase::PanoramaOptions opts;
00600 
00601         if (param.fisheye)
00602         {
00603             opts.setProjection(HuginBase::PanoramaOptions::FULL_FRAME_FISHEYE);
00604         }
00605         else
00606         {
00607             opts.setProjection(HuginBase::PanoramaOptions::RECTILINEAR);
00608         }
00609         opts.setHFOV(srcImg.getHFOV(), false);
00610         opts.setWidth(srcImg.getSize().x, false);
00611         opts.setHeight(srcImg.getSize().y);
00612         // output to tiff format
00613         opts.outputFormat = HuginBase::PanoramaOptions::TIFF_m;
00614         opts.tiff_saveROI = false;
00615         // m estimator, to be more robust against points on moving objects
00616         opts.huberSigma = 2;
00617         // save also exposure value of first image
00618         opts.outputExposureValue = srcImg.getExposureValue();
00619         pano.setOptions(opts);
00620 
00621         // variables that should be optimized
00622         // optimize nothing in the first image
00623         HuginBase::OptimizeVector optvars(1);
00624         std::vector<unsigned int> images;
00625         images.push_back(0);
00626 
00627         HuginBase::StandardImageVariableGroups variable_groups(pano);
00628 
00629         // loop to add images.
00630         for (int i = 1; i < (int) files.size(); i++)
00631         {
00632             // add next image.
00633             srcImg.setFilename(files[i]);
00634             srcImg.readEXIF();
00635             srcImg.applyEXIFValues();
00636             if (srcImg.getSize().x == 0 || srcImg.getSize().y == 0)
00637             {
00638                 std::cerr << "Could not decode image: " << files[i] << "Unsupported image file format" << std::endl;
00639                 return 1;
00640             }
00641             if (pano.getImage(0).getSize() != srcImg.getSize())
00642             {
00643                 std::cerr << "Images have different sizes." << std::endl
00644                     << files[0] << " has " << pano.getImage(0).getSize() << " pixel, while " << std::endl
00645                     << files[i] << " has " << srcImg.getSize() << " pixel." << std::endl
00646                     << "This is not supported. Align_image_stack works only with images of the same size." << std::endl;
00647                 return 1;
00648             };
00649 
00650             if(param.loadDistortion)
00651             {
00652                 if(srcImg.readDistortionFromDB())
00653                 {
00654                     std::cout << "\tRead distortion data from lens database." << std::endl;
00655                 }
00656                 else
00657                 {
00658                     std::cout << "\tNo valid distortion data found in lens database." << std::endl;
00659                 }
00660             }
00661 
00662             if (param.hfov > 0)
00663             {
00664                 srcImg.setHFOV(param.hfov);
00665             }
00666             else if (srcImg.getCropFactor() == 0)
00667             {
00668                 // could not read HFOV, assuming default: 50
00669                 srcImg.setHFOV(50);
00670             }
00671 
00672             int imgNr = pano.addImage(srcImg);
00673             variable_groups.update();
00674             // each image shares the same lens.
00675             variable_groups.getLenses().switchParts(imgNr, 0);
00676             // unlink HFOV?
00677             if (param.optHFOV)
00678             {
00679                 pano.unlinkImageVariableHFOV(0);
00680             }
00681             if (param.optDistortion)
00682             {
00683                 pano.unlinkImageVariableRadialDistortion(0);
00684             }
00685             if (param.optCenter)
00686             {
00687                 pano.unlinkImageVariableRadialDistortionCenterShift(0);
00688             }
00689             // All images are in the same stack: Link the stack variable.
00690             pano.linkImageVariableStack(imgNr, 0);
00691             // exposure value is linked, reset to value found in EXIF
00692             pano.unlinkImageVariableExposureValue(0);
00693             srcImg.setExposureValue(srcImg.calcExifExposureValue());
00694             pano.setSrcImage(i, srcImg);
00695             images.push_back(i);
00696             // optimize yaw, roll, pitch
00697             std::set<std::string> vars;
00698             vars.insert("y");
00699             vars.insert("p");
00700             vars.insert("r");
00701             if (param.optHFOV)
00702             {
00703                 vars.insert("v");
00704             }
00705             if (param.optDistortion)
00706             {
00707                 vars.insert("a");
00708                 vars.insert("b");
00709                 vars.insert("c");
00710             }
00711             if (param.optCenter)
00712             {
00713                 vars.insert("d");
00714                 vars.insert("e");
00715             }
00716             if (param.optX)
00717             {
00718                 vars.insert("TrX");
00719             }
00720             if (param.optY)
00721             {
00722                 vars.insert("TrY");
00723             }
00724             if (param.optZ)
00725             {
00726                 vars.insert("TrZ");
00727             }
00728             optvars.push_back(vars);
00729 
00730         };
00731 
00732         // sort now all images by exposure value
00733         if (param.sortImagesByEv)
00734         {
00735             std::sort(images.begin(), images.end(), SortImageVectorEV(&pano));
00736         };
00737 
00738         // load first image
00739         vigra::ImageImportInfo firstImgInfo(pano.getSrcImage(images[0]).getFilename().c_str());
00740 
00741         // original size
00742         ImageType* leftImgOrig = new ImageType(firstImgInfo.size());
00743         // rescale image
00744         ImageType* leftImg = new ImageType();
00745         {
00746             if(firstImgInfo.numExtraBands() == 1)
00747             {
00748                 vigra::BImage alpha(firstImgInfo.size());
00749                 vigra::importImageAlpha(firstImgInfo, destImage(*leftImgOrig), destImage(alpha));
00750             }
00751             else if (firstImgInfo.numExtraBands() == 0)
00752             {
00753                 vigra::importImage(firstImgInfo, destImage(*leftImgOrig));
00754             }
00755             else
00756             {
00757                 vigra_fail("Images with multiple extra (alpha) channels not supported");
00758             }
00759             vigra_ext::reduceNTimes(*leftImgOrig, *leftImg, param.pyrLevel);
00760         }
00761 
00762 
00763         ImageType* rightImg = new ImageType(leftImg->size());
00764         ImageType* rightImgOrig = new ImageType(leftImgOrig->size());
00765         // loop to add control points between them.
00766         for (int i = 1; i < (int) images.size(); i++)
00767         {
00768             if (g_verbose > 0)
00769             {
00770                 std::cout << "Creating control points between " << pano.getSrcImage(images[i-1]).getFilename().c_str() << " and " <<
00771                      pano.getSrcImage(images[i]).getFilename().c_str() << std::endl;
00772             }
00773 
00774             // load the actual image data of the next image
00775             vigra::ImageImportInfo nextImgInfo(pano.getSrcImage(images[i]).getFilename().c_str());
00776             assert(nextImgInfo.size() == firstImgInfo.size());
00777             {
00778                 if (nextImgInfo.numExtraBands() == 1)
00779                 {
00780                     vigra::BImage alpha(nextImgInfo.size());
00781                     vigra::importImageAlpha(nextImgInfo, destImage(*rightImgOrig), destImage(alpha));
00782                 }
00783                 else if (nextImgInfo.numExtraBands() == 0)
00784                 {
00785                     vigra::importImage(nextImgInfo, destImage(*rightImgOrig));
00786                 }
00787                 else
00788                 {
00789                     vigra_fail("Images with multiple extra (alpha) channels not supported");
00790                 }
00791                 vigra_ext::reduceNTimes(*rightImgOrig, *rightImg, param.pyrLevel);
00792             }
00793 
00794             // add control points.
00795             // work on smaller images
00796             // TODO: or use a fast interest point operator.
00797             createCtrlPoints(pano, images[i-1], *leftImg, *leftImgOrig, images[i], *rightImg, *rightImgOrig, param.pyrLevel, 2, param.nPoints, param.grid, param.corrThresh, param.stereo);
00798 
00799             // swap images;
00800             delete leftImg;
00801             delete leftImgOrig;
00802             leftImg = rightImg;
00803             leftImgOrig = rightImgOrig;
00804             rightImg = new ImageType(leftImg->size());
00805             rightImgOrig = new ImageType(leftImgOrig->size());
00806         }
00807         delete leftImg;
00808         delete rightImg;
00809         delete leftImgOrig;
00810         delete rightImgOrig;
00811 
00812         // optimize everything.
00813         pano.setOptimizeVector(optvars);
00814         // disable optimizer progress messages if -v not given
00815         if (g_verbose == 0)
00816         {
00817             PT_setProgressFcn(ptProgress);
00818             PT_setInfoDlgFcn(ptinfoDlg);
00819         };
00820         bool optimizeError = (HuginBase::PTools::optimize(pano) > 0);
00821 
00822         // need to do some basic outlier pruning.
00823         // remove all points with error higher than a specified threshold
00824         if (param.cpErrorThreshold > 0)
00825         {
00826             HuginBase::CPVector cps = pano.getCtrlPoints();
00827             HuginBase::CPVector newCPs;
00828             for (int i=0; i < (int)cps.size(); i++)
00829             {
00830                 if (cps[i].error < param.cpErrorThreshold ||
00831                     cps[i].mode == HuginBase::ControlPoint::X)   // preserve the vertical control point for stereo alignment
00832                 {
00833                     newCPs.push_back(cps[i]);
00834                 }
00835             }
00836             if (g_verbose > 0)
00837             {
00838                 std::cout << "Ctrl points before pruning: " << cps.size() << ", after: " << newCPs.size() << std::endl;
00839             }
00840             pano.setCtrlPoints(newCPs);
00841             if (param.stereo_window)
00842             {
00843                 alignStereoWindow(pano, param.pop_out);
00844             }
00845             // reoptimize
00846             optimizeError = (HuginBase::PTools::optimize(pano) > 0);
00847         }
00848 
00849         if (param.crop)
00850         {
00851             autoCrop(pano);
00852         }
00853 
00854         HuginBase::UIntSet imgs = pano.getActiveImages();
00855         if (optimizeError)
00856         {
00857             if (param.ptoFile.size() > 0)
00858             {
00859                 std::ofstream script(param.ptoFile.c_str());
00860                 pano.printPanoramaScript(script, optvars, pano.getOptions(), imgs, false, "");
00861             }
00862             std::cerr << "An error occurred during optimization." << std::endl;
00863             std::cerr << "Try adding \"-p debug.pto\" and checking output." << std::endl;
00864             std::cerr << "Exiting..." << std::endl;
00865             return 1;
00866         }
00867 
00868         if (param.hdrFile.size())
00869         {
00870             // TODO: photometric alignment (HDR, fixed white balance)
00871             //utils::StreamProgressReporter progress(2.0);
00872             //loadImgsAndExtractPoints(pano, nPoints, pyrLevel, randomPoints, progress, points);
00873             //smartOptimizePhotometric
00874 
00875             // switch to HDR output mode
00876             HuginBase::PanoramaOptions opts = pano.getOptions();
00877             opts.outputFormat = HuginBase::PanoramaOptions::HDR;
00878             opts.outputPixelType = "FLOAT";
00879             opts.outputMode = HuginBase::PanoramaOptions::OUTPUT_HDR;
00880             opts.outfile = param.hdrFile;
00881             pano.setOptions(opts);
00882 
00883             // remap all images
00884             AppBase::ProgressDisplay* progress;
00885             if(g_verbose > 0)
00886             {
00887                 progress = new AppBase::StreamProgressDisplay(std::cout);
00888             }
00889             else
00890             {
00891                 progress = new AppBase::DummyProgressDisplay();
00892             };
00893             HuginBase::Nona::stitchPanorama(pano, pano.getOptions(),
00894                            progress, opts.outfile, imgs);
00895             std::cout << "Written HDR output to " << opts.outfile << std::endl;
00896             delete progress;
00897         }
00898         if (param.alignedPrefix.size())
00899         {
00900             // disable all exposure compensation stuff.
00901             HuginBase::PanoramaOptions opts = pano.getOptions();
00902             opts.outputMode = HuginBase::PanoramaOptions::OUTPUT_LDR;
00903             opts.outputFormat = HuginBase::PanoramaOptions::TIFF_m;
00904             opts.outputPixelType = "";
00905             opts.outfile = param.alignedPrefix;
00906             opts.remapUsingGPU = param.gpu;
00907             pano.setOptions(opts);
00908             // remap all images
00909             AppBase::ProgressDisplay* progress;
00910             if(g_verbose > 0)
00911             {
00912                 progress = new AppBase::StreamProgressDisplay(std::cout);
00913             }
00914             else
00915             {
00916                 progress = new AppBase::DummyProgressDisplay();
00917             };
00918             // pass option to ignore exposure to stitcher
00919             HuginBase::Nona::AdvancedOptions advOpts;
00920             HuginBase::Nona::SetAdvancedOption(advOpts, "ignoreExposure", true);
00921             HuginBase::Nona::stitchPanorama(pano, pano.getOptions(),
00922                            progress, opts.outfile, imgs, advOpts);
00923             delete progress;
00924             std::cout << "Written aligned images to files with prefix \"" << opts.outfile << "\"" << std::endl;
00925         }
00926 
00927         // At this point we have panorama options set according to the output
00928         if (param.ptoFile.size() > 0)
00929         {
00930             std::ofstream script(param.ptoFile.c_str());
00931             pano.printPanoramaScript(script, optvars, pano.getOptions(), imgs, false, "");
00932             std::cout << "Written project file " << param.ptoFile << std::endl;
00933         }
00934 
00935 
00936     }
00937     catch (std::exception& e)
00938     {
00939         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
00940         return 1;
00941     }
00942     return 0;
00943 }
00944 
00945 int main(int argc, char* argv[])
00946 {
00947     // parse arguments
00948     const char* optstring = "a:ef:g:hlmdiSAPCp:vo:s:t:c:xyz";
00949     int c;
00950 
00951     g_verbose = 0;
00952 
00953     Parameters param;
00954     enum
00955     {
00956         CORRTHRESH=1000,
00957         THREADS,
00958         GPU,
00959         LENSDB,
00960         USEGIVENORDER,
00961     };
00962 
00963     static struct option longOptions[] =
00964     {
00965         {"corr", required_argument, NULL, CORRTHRESH },
00966         {"threads", required_argument, NULL, THREADS },
00967         {"gpu", no_argument, NULL, GPU },
00968         {"distortion", no_argument, NULL, LENSDB },
00969         {"use-given-order", no_argument, NULL, USEGIVENORDER },
00970         {"help", no_argument, NULL, 'h' },
00971         0
00972     };
00973     while ((c = getopt_long(argc, argv, optstring, longOptions, nullptr)) != -1)
00974     {
00975         switch (c)
00976         {
00977             case 'a':
00978                 param.alignedPrefix = optarg;
00979                 break;
00980             case 'c':
00981                 param.nPoints = atoi(optarg);
00982                 if (param.nPoints < 1)
00983                 {
00984                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: Number of points/grid (-c) must be at least 1" << std::endl;
00985                     return 1;
00986                 }
00987                 break;
00988             case 'e':
00989                 param.fisheye = true;
00990                 break;
00991             case 'f':
00992                 param.hfov = atof(optarg);
00993                 if (param.hfov <= 0)
00994                 {
00995                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: HFOV (-f) must be greater than 0" << std::endl;
00996                     return 1;
00997                 }
00998                 break;
00999             case 'g':
01000                 param.grid = atoi(optarg);
01001                 if (param.grid < 1 || param.grid>50)
01002                 {
01003                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: number of grid cells (-g) should be between 1 and 50" << std::endl;
01004                     return 1;
01005                 }
01006                 break;
01007             case 'l':
01008                 param.linear = true;
01009                 break;
01010             case 'm':
01011                 param.optHFOV = true;
01012                 break;
01013             case 'd':
01014                 param.optDistortion = true;
01015                 break;
01016             case 'i':
01017                 param.optCenter = true;
01018                 break;
01019             case 'x':
01020                 param.optX = true;
01021                 break;
01022             case 'y':
01023                 param.optY = true;
01024                 break;
01025             case 'z':
01026                 param.optZ = true;
01027                 break;
01028             case 'S':
01029                 param.stereo = true;
01030                 break;
01031             case 'A':
01032                 param.stereo = true;
01033                 param.stereo_window = true;
01034                 break;
01035             case 'P':
01036                 param.stereo = true;
01037                 param.stereo_window = true;
01038                 param.pop_out = true;
01039                 break;
01040             case 'C':
01041                 param.crop = true;
01042                 break;
01043             case 't':
01044                 param.cpErrorThreshold = atof(optarg);
01045                 if (param.cpErrorThreshold <= 0)
01046                 {
01047                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: control point error threshold (-t) must be greater than 0" << std::endl;
01048                     return 1;
01049                 }
01050                 break;
01051             case 'p':
01052                 param.ptoFile = optarg;
01053                 break;
01054             case 'o':
01055                 param.hdrFile = optarg;
01056                 break;
01057             case 'v':
01058                 g_verbose++;
01059                 break;
01060             case 'h':
01061                 usage(hugin_utils::stripPath(argv[0]).c_str());
01062                 return 0;
01063             case 's':
01064                 param.pyrLevel = atoi(optarg);
01065                 if (param.pyrLevel < 0 || param.pyrLevel >8)
01066                 {
01067                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: scaling (-s) should be between 0 and 8" << std::endl;
01068                     return 1;
01069                 }
01070                 break;
01071             case CORRTHRESH:
01072                 param.corrThresh = atof(optarg);
01073                 if (param.corrThresh <= 0 || param.corrThresh > 1.0)
01074                 {
01075                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: correlation should be between 0 and 1" << endl;
01076                     return 1;
01077                 };
01078                 break;
01079             case THREADS:
01080                 std::cout << "WARNING: Switch --threads is deprecated. Set environment variable OMP_NUM_THREADS instead" << std::endl;
01081                 break;
01082             case GPU:
01083                 param.gpu = true;
01084                 break;
01085             case LENSDB:
01086                 param.loadDistortion = true;
01087                 break;
01088             case USEGIVENORDER:
01089                 param.sortImagesByEv = false;
01090                 break;
01091             case ':':
01092             case '?':
01093                 // missing argument or invalid switch
01094                 return 1;
01095                 break;
01096             default:
01097                 // this should not happen
01098                 abort();
01099         }
01100     }
01101 
01102     // use always given image order for stereo options
01103     if (param.stereo)
01104     {
01105         param.sortImagesByEv = false;
01106     };
01107     unsigned nFiles = argc - optind;
01108     if (nFiles < 2)
01109     {
01110         std::cerr << hugin_utils::stripPath(argv[0]) << ": At least two files need to be specified" << std::endl;
01111         return 1;
01112     }
01113 
01114     if (param.hdrFile.empty() && param.ptoFile.empty() && param.alignedPrefix.empty())
01115     {
01116         std::cerr << hugin_utils::stripPath(argv[0]) << ": Please specify at least one of the -p, -o or -a options." << std::endl;
01117         return 1;
01118     }
01119 
01120     // extract file names
01121     std::vector<std::string> files;
01122     for (size_t i=0; i < nFiles; i++)
01123     {
01124         files.push_back(argv[optind+i]);
01125     }
01126 
01127     std::string pixelType;
01128 
01129     bool grayscale = false;
01130     int returnValue = 1;
01131     try
01132     {
01133         vigra::ImageImportInfo firstImgInfo(files[0].c_str());
01134         pixelType = firstImgInfo.getPixelType();
01135         if (firstImgInfo.numExtraBands()>1)
01136         {
01137             std::cerr << "ERROR: images with several alpha channels are not supported." << std::endl;
01138             return 1;
01139         };
01140         grayscale = firstImgInfo.isGrayscale();
01141     }
01142     catch (std::exception& e)
01143     {
01144         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
01145         return 1;
01146     }
01147 
01148     if(param.gpu)
01149     {
01150         param.gpu=hugin_utils::initGPU(&argc, argv);
01151     };
01152     if (grayscale)
01153     {
01154         if (pixelType == "UINT8")
01155         {
01156             returnValue=main2<vigra::UInt8>(files, param);
01157         }
01158         else if (pixelType == "INT16")
01159         {
01160             returnValue=main2<vigra::Int16>(files, param);
01161         }
01162         else if (pixelType == "UINT16")
01163         {
01164             returnValue = main2<vigra::UInt16>(files, param);
01165         }
01166         else if (pixelType == "FLOAT")
01167         {
01168             returnValue=main2<float>(files, param);
01169         }
01170         else
01171         {
01172             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01173         }
01174     }
01175     else
01176     {
01177         if (pixelType == "UINT8")
01178         {
01179             returnValue = main2<vigra::RGBValue<vigra::UInt8> >(files, param);
01180         }
01181         else if (pixelType == "INT16")
01182         {
01183             returnValue = main2<vigra::RGBValue<vigra::Int16> >(files, param);
01184         }
01185         else if (pixelType == "UINT16")
01186         {
01187             returnValue = main2<vigra::RGBValue<vigra::UInt16> >(files, param);
01188         }
01189         else if (pixelType == "FLOAT")
01190         {
01191             returnValue = main2<vigra::RGBValue<float> >(files, param);
01192         }
01193         else
01194         {
01195             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01196         }
01197     };
01198 
01199     if(param.gpu)
01200     {
01201         hugin_utils::wrapupGPU();
01202     };
01203     HuginBase::LensDB::LensDB::Clean();
01204     return returnValue;
01205 }

Generated on 8 Dec 2016 for Hugintrunk by  doxygen 1.4.7