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

Generated on 25 Aug 2016 for Hugintrunk by  doxygen 1.4.7