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::cout << 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             // reset size to force new detection of image size
00635             srcImg.setSize(vigra::Size2D());
00636             srcImg.readEXIF();
00637             srcImg.applyEXIFValues();
00638             if (srcImg.getSize().x == 0 || srcImg.getSize().y == 0)
00639             {
00640                 std::cerr << "Could not decode image: " << files[i] << "Unsupported image file format" << std::endl;
00641                 return 1;
00642             }
00643             if (pano.getImage(0).getSize() != srcImg.getSize())
00644             {
00645                 std::cerr << "Images have different sizes." << std::endl
00646                     << files[0] << " has " << pano.getImage(0).getSize() << " pixel, while " << std::endl
00647                     << files[i] << " has " << srcImg.getSize() << " pixel." << std::endl
00648                     << "This is not supported. Align_image_stack works only with images of the same size." << std::endl;
00649                 return 1;
00650             };
00651 
00652             if(param.loadDistortion)
00653             {
00654                 if(srcImg.readDistortionFromDB())
00655                 {
00656                     std::cout << "\tRead distortion data from lens database." << std::endl;
00657                 }
00658                 else
00659                 {
00660                     std::cout << "\tNo valid distortion data found in lens database." << std::endl;
00661                 }
00662             }
00663 
00664             if (param.hfov > 0)
00665             {
00666                 srcImg.setHFOV(param.hfov);
00667             }
00668             else if (srcImg.getCropFactor() == 0)
00669             {
00670                 // could not read HFOV, assuming default: 50
00671                 srcImg.setHFOV(50);
00672             }
00673 
00674             int imgNr = pano.addImage(srcImg);
00675             variable_groups.update();
00676             // each image shares the same lens.
00677             variable_groups.getLenses().switchParts(imgNr, 0);
00678             // unlink HFOV?
00679             if (param.optHFOV)
00680             {
00681                 pano.unlinkImageVariableHFOV(0);
00682             }
00683             if (param.optDistortion)
00684             {
00685                 pano.unlinkImageVariableRadialDistortion(0);
00686             }
00687             if (param.optCenter)
00688             {
00689                 pano.unlinkImageVariableRadialDistortionCenterShift(0);
00690             }
00691             // All images are in the same stack: Link the stack variable.
00692             pano.linkImageVariableStack(imgNr, 0);
00693             // exposure value is linked, reset to value found in EXIF
00694             pano.unlinkImageVariableExposureValue(0);
00695             srcImg.setExposureValue(srcImg.calcExifExposureValue());
00696             pano.setSrcImage(i, srcImg);
00697             images.push_back(i);
00698             // optimize yaw, roll, pitch
00699             std::set<std::string> vars;
00700             vars.insert("y");
00701             vars.insert("p");
00702             vars.insert("r");
00703             if (param.optHFOV)
00704             {
00705                 vars.insert("v");
00706             }
00707             if (param.optDistortion)
00708             {
00709                 vars.insert("a");
00710                 vars.insert("b");
00711                 vars.insert("c");
00712             }
00713             if (param.optCenter)
00714             {
00715                 vars.insert("d");
00716                 vars.insert("e");
00717             }
00718             if (param.optX)
00719             {
00720                 vars.insert("TrX");
00721             }
00722             if (param.optY)
00723             {
00724                 vars.insert("TrY");
00725             }
00726             if (param.optZ)
00727             {
00728                 vars.insert("TrZ");
00729             }
00730             optvars.push_back(vars);
00731 
00732         };
00733 
00734         // sort now all images by exposure value
00735         if (param.sortImagesByEv)
00736         {
00737             std::sort(images.begin(), images.end(), SortImageVectorEV(&pano));
00738         };
00739 
00740         // load first image
00741         vigra::ImageImportInfo firstImgInfo(pano.getSrcImage(images[0]).getFilename().c_str());
00742 
00743         // original size
00744         ImageType* leftImgOrig = new ImageType(firstImgInfo.size());
00745         // rescale image
00746         ImageType* leftImg = new ImageType();
00747         {
00748             if(firstImgInfo.numExtraBands() == 1)
00749             {
00750                 vigra::BImage alpha(firstImgInfo.size());
00751                 vigra::importImageAlpha(firstImgInfo, destImage(*leftImgOrig), destImage(alpha));
00752             }
00753             else if (firstImgInfo.numExtraBands() == 0)
00754             {
00755                 vigra::importImage(firstImgInfo, destImage(*leftImgOrig));
00756             }
00757             else
00758             {
00759                 vigra_fail("Images with multiple extra (alpha) channels not supported");
00760             }
00761             vigra_ext::reduceNTimes(*leftImgOrig, *leftImg, param.pyrLevel);
00762         }
00763 
00764 
00765         ImageType* rightImg = new ImageType(leftImg->size());
00766         ImageType* rightImgOrig = new ImageType(leftImgOrig->size());
00767         // loop to add control points between them.
00768         for (int i = 1; i < (int) images.size(); i++)
00769         {
00770             if (g_verbose > 0)
00771             {
00772                 std::cout << "Creating control points between " << pano.getSrcImage(images[i-1]).getFilename().c_str() << " and " <<
00773                      pano.getSrcImage(images[i]).getFilename().c_str() << std::endl;
00774             }
00775 
00776             // load the actual image data of the next image
00777             vigra::ImageImportInfo nextImgInfo(pano.getSrcImage(images[i]).getFilename().c_str());
00778             assert(nextImgInfo.size() == firstImgInfo.size());
00779             {
00780                 if (nextImgInfo.numExtraBands() == 1)
00781                 {
00782                     vigra::BImage alpha(nextImgInfo.size());
00783                     vigra::importImageAlpha(nextImgInfo, destImage(*rightImgOrig), destImage(alpha));
00784                 }
00785                 else if (nextImgInfo.numExtraBands() == 0)
00786                 {
00787                     vigra::importImage(nextImgInfo, destImage(*rightImgOrig));
00788                 }
00789                 else
00790                 {
00791                     vigra_fail("Images with multiple extra (alpha) channels not supported");
00792                 }
00793                 vigra_ext::reduceNTimes(*rightImgOrig, *rightImg, param.pyrLevel);
00794             }
00795 
00796             // add control points.
00797             // work on smaller images
00798             // TODO: or use a fast interest point operator.
00799             createCtrlPoints(pano, images[i-1], *leftImg, *leftImgOrig, images[i], *rightImg, *rightImgOrig, param.pyrLevel, 2, param.nPoints, param.grid, param.corrThresh, param.stereo);
00800 
00801             // swap images;
00802             delete leftImg;
00803             delete leftImgOrig;
00804             leftImg = rightImg;
00805             leftImgOrig = rightImgOrig;
00806             rightImg = new ImageType(leftImg->size());
00807             rightImgOrig = new ImageType(leftImgOrig->size());
00808         }
00809         delete leftImg;
00810         delete rightImg;
00811         delete leftImgOrig;
00812         delete rightImgOrig;
00813 
00814         // optimize everything.
00815         pano.setOptimizeVector(optvars);
00816         // disable optimizer progress messages if -v not given
00817         if (g_verbose == 0)
00818         {
00819             PT_setProgressFcn(ptProgress);
00820             PT_setInfoDlgFcn(ptinfoDlg);
00821         };
00822         bool optimizeError = (HuginBase::PTools::optimize(pano) > 0);
00823 
00824         // need to do some basic outlier pruning.
00825         // remove all points with error higher than a specified threshold
00826         if (param.cpErrorThreshold > 0)
00827         {
00828             HuginBase::CPVector cps = pano.getCtrlPoints();
00829             HuginBase::CPVector newCPs;
00830             for (int i=0; i < (int)cps.size(); i++)
00831             {
00832                 if (cps[i].error < param.cpErrorThreshold ||
00833                     cps[i].mode == HuginBase::ControlPoint::X)   // preserve the vertical control point for stereo alignment
00834                 {
00835                     newCPs.push_back(cps[i]);
00836                 }
00837             }
00838             if (g_verbose > 0)
00839             {
00840                 std::cout << "Ctrl points before pruning: " << cps.size() << ", after: " << newCPs.size() << std::endl;
00841             }
00842             pano.setCtrlPoints(newCPs);
00843             if (param.stereo_window)
00844             {
00845                 alignStereoWindow(pano, param.pop_out);
00846             }
00847             // reoptimize
00848             optimizeError = (HuginBase::PTools::optimize(pano) > 0);
00849         }
00850 
00851         if (param.crop)
00852         {
00853             autoCrop(pano);
00854         }
00855 
00856         HuginBase::UIntSet imgs = pano.getActiveImages();
00857         if (optimizeError)
00858         {
00859             if (param.ptoFile.size() > 0)
00860             {
00861                 std::ofstream script(param.ptoFile.c_str());
00862                 pano.printPanoramaScript(script, optvars, pano.getOptions(), imgs, false, "");
00863             }
00864             std::cerr << "An error occurred during optimization." << std::endl;
00865             std::cerr << "Try adding \"-p debug.pto\" and checking output." << std::endl;
00866             std::cerr << "Exiting..." << std::endl;
00867             return 1;
00868         }
00869 
00870         if (param.hdrFile.size())
00871         {
00872             // TODO: photometric alignment (HDR, fixed white balance)
00873             //utils::StreamProgressReporter progress(2.0);
00874             //loadImgsAndExtractPoints(pano, nPoints, pyrLevel, randomPoints, progress, points);
00875             //smartOptimizePhotometric
00876 
00877             // switch to HDR output mode
00878             HuginBase::PanoramaOptions opts = pano.getOptions();
00879             opts.outputFormat = HuginBase::PanoramaOptions::HDR;
00880             opts.outputPixelType = "FLOAT";
00881             opts.outputMode = HuginBase::PanoramaOptions::OUTPUT_HDR;
00882             opts.outfile = param.hdrFile;
00883             pano.setOptions(opts);
00884 
00885             // remap all images
00886             AppBase::ProgressDisplay* progress;
00887             if(g_verbose > 0)
00888             {
00889                 progress = new AppBase::StreamProgressDisplay(std::cout);
00890             }
00891             else
00892             {
00893                 progress = new AppBase::DummyProgressDisplay();
00894             };
00895             HuginBase::Nona::stitchPanorama(pano, pano.getOptions(),
00896                            progress, opts.outfile, imgs);
00897             std::cout << "Written HDR output to " << opts.outfile << std::endl;
00898             delete progress;
00899         }
00900         if (param.alignedPrefix.size())
00901         {
00902             // disable all exposure compensation stuff.
00903             HuginBase::PanoramaOptions opts = pano.getOptions();
00904             opts.outputMode = HuginBase::PanoramaOptions::OUTPUT_LDR;
00905             opts.outputFormat = HuginBase::PanoramaOptions::TIFF_m;
00906             opts.outputPixelType = "";
00907             opts.outfile = param.alignedPrefix;
00908             opts.remapUsingGPU = param.gpu;
00909             pano.setOptions(opts);
00910             // remap all images
00911             AppBase::ProgressDisplay* progress;
00912             if(g_verbose > 0)
00913             {
00914                 progress = new AppBase::StreamProgressDisplay(std::cout);
00915             }
00916             else
00917             {
00918                 progress = new AppBase::DummyProgressDisplay();
00919             };
00920             // pass option to ignore exposure to stitcher
00921             HuginBase::Nona::AdvancedOptions advOpts;
00922             HuginBase::Nona::SetAdvancedOption(advOpts, "ignoreExposure", true);
00923             HuginBase::Nona::stitchPanorama(pano, pano.getOptions(),
00924                            progress, opts.outfile, imgs, advOpts);
00925             delete progress;
00926             std::cout << "Written aligned images to files with prefix \"" << opts.outfile << "\"" << std::endl;
00927         }
00928 
00929         // At this point we have panorama options set according to the output
00930         if (param.ptoFile.size() > 0)
00931         {
00932             std::ofstream script(param.ptoFile.c_str());
00933             pano.printPanoramaScript(script, optvars, pano.getOptions(), imgs, false, "");
00934             std::cout << "Written project file " << param.ptoFile << std::endl;
00935         }
00936 
00937 
00938     }
00939     catch (std::exception& e)
00940     {
00941         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
00942         return 1;
00943     }
00944     return 0;
00945 }
00946 
00947 int main(int argc, char* argv[])
00948 {
00949     // parse arguments
00950     const char* optstring = "a:ef:g:hlmdiSAPCp:vo:s:t:c:xyz";
00951     int c;
00952 
00953     g_verbose = 0;
00954 
00955     Parameters param;
00956     enum
00957     {
00958         CORRTHRESH=1000,
00959         THREADS,
00960         GPU,
00961         LENSDB,
00962         USEGIVENORDER,
00963     };
00964 
00965     static struct option longOptions[] =
00966     {
00967         {"corr", required_argument, NULL, CORRTHRESH },
00968         {"threads", required_argument, NULL, THREADS },
00969         {"gpu", no_argument, NULL, GPU },
00970         {"distortion", no_argument, NULL, LENSDB },
00971         {"use-given-order", no_argument, NULL, USEGIVENORDER },
00972         {"help", no_argument, NULL, 'h' },
00973         0
00974     };
00975     while ((c = getopt_long(argc, argv, optstring, longOptions, nullptr)) != -1)
00976     {
00977         switch (c)
00978         {
00979             case 'a':
00980                 param.alignedPrefix = optarg;
00981                 break;
00982             case 'c':
00983                 param.nPoints = atoi(optarg);
00984                 if (param.nPoints < 1)
00985                 {
00986                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: Number of points/grid (-c) must be at least 1" << std::endl;
00987                     return 1;
00988                 }
00989                 break;
00990             case 'e':
00991                 param.fisheye = true;
00992                 break;
00993             case 'f':
00994                 param.hfov = atof(optarg);
00995                 if (param.hfov <= 0)
00996                 {
00997                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: HFOV (-f) must be greater than 0" << std::endl;
00998                     return 1;
00999                 }
01000                 break;
01001             case 'g':
01002                 param.grid = atoi(optarg);
01003                 if (param.grid < 1 || param.grid>50)
01004                 {
01005                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: number of grid cells (-g) should be between 1 and 50" << std::endl;
01006                     return 1;
01007                 }
01008                 break;
01009             case 'l':
01010                 param.linear = true;
01011                 break;
01012             case 'm':
01013                 param.optHFOV = true;
01014                 break;
01015             case 'd':
01016                 param.optDistortion = true;
01017                 break;
01018             case 'i':
01019                 param.optCenter = true;
01020                 break;
01021             case 'x':
01022                 param.optX = true;
01023                 break;
01024             case 'y':
01025                 param.optY = true;
01026                 break;
01027             case 'z':
01028                 param.optZ = true;
01029                 break;
01030             case 'S':
01031                 param.stereo = true;
01032                 break;
01033             case 'A':
01034                 param.stereo = true;
01035                 param.stereo_window = true;
01036                 break;
01037             case 'P':
01038                 param.stereo = true;
01039                 param.stereo_window = true;
01040                 param.pop_out = true;
01041                 break;
01042             case 'C':
01043                 param.crop = true;
01044                 break;
01045             case 't':
01046                 param.cpErrorThreshold = atof(optarg);
01047                 if (param.cpErrorThreshold <= 0)
01048                 {
01049                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: control point error threshold (-t) must be greater than 0" << std::endl;
01050                     return 1;
01051                 }
01052                 break;
01053             case 'p':
01054                 param.ptoFile = optarg;
01055                 break;
01056             case 'o':
01057                 param.hdrFile = optarg;
01058                 break;
01059             case 'v':
01060                 g_verbose++;
01061                 break;
01062             case 'h':
01063                 usage(hugin_utils::stripPath(argv[0]).c_str());
01064                 return 0;
01065             case 's':
01066                 param.pyrLevel = atoi(optarg);
01067                 if (param.pyrLevel < 0 || param.pyrLevel >8)
01068                 {
01069                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: scaling (-s) should be between 0 and 8" << std::endl;
01070                     return 1;
01071                 }
01072                 break;
01073             case CORRTHRESH:
01074                 param.corrThresh = atof(optarg);
01075                 if (param.corrThresh <= 0 || param.corrThresh > 1.0)
01076                 {
01077                     std::cerr << hugin_utils::stripPath(argv[0]) << ": Invalid parameter: correlation should be between 0 and 1" << endl;
01078                     return 1;
01079                 };
01080                 break;
01081             case THREADS:
01082                 std::cout << "WARNING: Switch --threads is deprecated. Set environment variable OMP_NUM_THREADS instead" << std::endl;
01083                 break;
01084             case GPU:
01085                 param.gpu = true;
01086                 break;
01087             case LENSDB:
01088                 param.loadDistortion = true;
01089                 break;
01090             case USEGIVENORDER:
01091                 param.sortImagesByEv = false;
01092                 break;
01093             case ':':
01094             case '?':
01095                 // missing argument or invalid switch
01096                 return 1;
01097                 break;
01098             default:
01099                 // this should not happen
01100                 abort();
01101         }
01102     }
01103 
01104     // use always given image order for stereo options
01105     if (param.stereo)
01106     {
01107         param.sortImagesByEv = false;
01108     };
01109     unsigned nFiles = argc - optind;
01110     if (nFiles < 2)
01111     {
01112         std::cerr << hugin_utils::stripPath(argv[0]) << ": At least two files need to be specified" << std::endl;
01113         return 1;
01114     }
01115 
01116     if (param.hdrFile.empty() && param.ptoFile.empty() && param.alignedPrefix.empty())
01117     {
01118         std::cerr << hugin_utils::stripPath(argv[0]) << ": Please specify at least one of the -p, -o or -a options." << std::endl;
01119         return 1;
01120     }
01121 
01122     // extract file names
01123     std::vector<std::string> files;
01124     for (size_t i=0; i < nFiles; i++)
01125     {
01126         files.push_back(argv[optind+i]);
01127     }
01128 
01129     std::string pixelType;
01130 
01131     bool grayscale = false;
01132     int returnValue = 1;
01133     try
01134     {
01135         vigra::ImageImportInfo firstImgInfo(files[0].c_str());
01136         pixelType = firstImgInfo.getPixelType();
01137         if (firstImgInfo.numExtraBands()>1)
01138         {
01139             std::cerr << "ERROR: images with several alpha channels are not supported." << std::endl;
01140             return 1;
01141         };
01142         grayscale = firstImgInfo.isGrayscale();
01143     }
01144     catch (std::exception& e)
01145     {
01146         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
01147         return 1;
01148     }
01149 
01150     if(param.gpu)
01151     {
01152         param.gpu=hugin_utils::initGPU(&argc, argv);
01153     };
01154     if (grayscale)
01155     {
01156         if (pixelType == "UINT8")
01157         {
01158             returnValue=main2<vigra::UInt8>(files, param);
01159         }
01160         else if (pixelType == "INT16")
01161         {
01162             returnValue=main2<vigra::Int16>(files, param);
01163         }
01164         else if (pixelType == "UINT16")
01165         {
01166             returnValue = main2<vigra::UInt16>(files, param);
01167         }
01168         else if (pixelType == "FLOAT")
01169         {
01170             returnValue=main2<float>(files, param);
01171         }
01172         else
01173         {
01174             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01175         }
01176     }
01177     else
01178     {
01179         if (pixelType == "UINT8")
01180         {
01181             returnValue = main2<vigra::RGBValue<vigra::UInt8> >(files, param);
01182         }
01183         else if (pixelType == "INT16")
01184         {
01185             returnValue = main2<vigra::RGBValue<vigra::Int16> >(files, param);
01186         }
01187         else if (pixelType == "UINT16")
01188         {
01189             returnValue = main2<vigra::RGBValue<vigra::UInt16> >(files, param);
01190         }
01191         else if (pixelType == "FLOAT")
01192         {
01193             returnValue = main2<vigra::RGBValue<float> >(files, param);
01194         }
01195         else
01196         {
01197             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01198         }
01199     };
01200 
01201     if(param.gpu)
01202     {
01203         hugin_utils::wrapupGPU();
01204     };
01205     HuginBase::LensDB::LensDB::Clean();
01206     return returnValue;
01207 }

Generated on 24 Nov 2017 for Hugintrunk by  doxygen 1.4.7