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

Generated on 31 Jul 2015 for Hugintrunk by  doxygen 1.4.7