tca_correct.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_ext/impexalpha.hxx>
00034 #include <vigra/cornerdetection.hxx>
00035 #include <vigra/localminmax.hxx>
00036 #include <hugin_utils/utils.h>
00037 #include <hugin_math/hugin_math.h>
00038 
00039 #include "vigra/stdimage.hxx"
00040 #include "vigra/stdimagefunctions.hxx"
00041 #include "vigra/functorexpression.hxx"
00042 #include "vigra/transformimage.hxx"
00043 
00044 #include <vigra_ext/Pyramid.h>
00045 #include <vigra_ext/Correlation.h>
00046 #include <vigra_ext/InterestPoints.h>
00047 #include <vigra_ext/utils.h>
00048 
00049 #include <panodata/Panorama.h>
00050 #include <panodata/StandardImageVariableGroups.h>
00051 #include <panotools/PanoToolsOptimizerWrapper.h>
00052 #include <algorithms/optimizer/PTOptimizer.h>
00053 #include <nona/Stitcher.h>
00054 #include <foreign/levmar/levmar.h>
00055 #include <hugin_utils/openmp_lock.h>
00056 #include <lensdb/LensDB.h>
00057 
00058 #include <getopt.h>
00059 #ifndef _WIN32
00060 #include <unistd.h>
00061 #endif
00062 
00063 #include <tiff.h>
00064 
00065 #define DEFAULT_OPTIMISATION_PARAMETER "abcvde"
00066 
00067 struct Parameters
00068 {
00069     Parameters()
00070     {
00071         cpErrorThreshold = 1.5;
00072         optMethod = 0;
00073         load = false;
00074         reset = false;
00075         saveDB = false;
00076         scale=2;
00077         nPoints=10;
00078         grid = 10;
00079         verbose = 0;
00080     }
00081 
00082     double cpErrorThreshold;
00083     int optMethod;
00084     bool load;
00085     bool reset;
00086     bool saveDB;
00087     std::set<std::string> optvars;
00088 
00089     std::string alignedPrefix;
00090     std::string ptoFile;
00091     std::string ptoOutputFile;
00092     std::string basename;
00093 
00094     std::string red_name;
00095     std::string green_name;
00096     std::string blue_name;
00097 
00098     double scale;
00099     int nPoints;
00100     unsigned grid;
00101     int verbose;
00102 };
00103 
00104 Parameters g_param;
00105 
00106 // Optimiser code
00107 struct OptimData
00108 {
00109     HuginBase::PanoramaData& m_pano;
00110     double huberSigma;
00111     const HuginBase::OptimizeVector& m_optvars;
00112 
00113     double m_dist[3][3]; // a,b,c for all imgs
00114     double m_shift[2];   // x,y shift
00115     double m_hfov[3];
00116     double m_center[2];  // center of image (without shift)
00117     std::vector<double*> m_mapping;
00118 
00119     int m_maxIter;
00120 
00121     OptimData(HuginBase::PanoramaData& pano, const HuginBase::OptimizeVector& optvars,
00122         double mEstimatorSigma, int maxIter)
00123         : m_pano(pano), huberSigma(mEstimatorSigma), m_optvars(optvars), m_maxIter(maxIter)
00124     {
00125         assert(m_pano.getNrOfImages() == m_optvars.size());
00126         assert(m_pano.getNrOfImages() == 3);
00127         LoadFromImgs();
00128 
00129         for (unsigned int i = 0; i<3; i++)
00130         {
00131             const std::set<std::string> vars = m_optvars[i];
00132             for (std::set<std::string>::const_iterator it = vars.begin(); it != vars.end(); ++it)
00133             {
00134                 const char var = (*it)[0];
00135                 if ((var >= 'a') && (var <= 'c'))
00136                 {
00137                     m_mapping.push_back(&(m_dist[i][var - 'a']));
00138                 }
00139                 else if ((var == 'd') || (var == 'e'))
00140                 {
00141                     m_mapping.push_back(&(m_shift[var - 'd']));
00142                 }
00143                 else if (var == 'v')
00144                 {
00145                     m_mapping.push_back(&(m_hfov[i]));
00146                 }
00147                 else
00148                 {
00149                     std::cerr << "Unknown parameter detected, ignoring!" << std::endl;
00150                 }
00151             }
00152         }
00153     }
00154 
00156     void ToX(double* x)
00157     {
00158         for (size_t i=0; i < m_mapping.size(); i++)
00159         {
00160             x[i] = *(m_mapping[i]);
00161         }
00162     }
00163 
00165     void FromX(double* x)
00166     {
00167         for (size_t i=0; i < m_mapping.size(); i++)
00168         {
00169             *(m_mapping[i]) = x[i];
00170         }
00171     }
00172 
00173     void LoadFromImgs()
00174     {
00175         for (unsigned int i = 0; i < 3; i++)
00176         {
00177             HuginBase::SrcPanoImage img = m_pano.getSrcImage(i);
00178             m_hfov[i] = img.getHFOV();
00179             m_dist[i][0] = img.getRadialDistortion()[0];
00180             m_dist[i][1] = img.getRadialDistortion()[1];
00181             m_dist[i][2] = img.getRadialDistortion()[2];
00182             if (i == 0)
00183             {
00184                 m_shift[0] = img.getRadialDistortionCenterShift().x;
00185                 m_shift[1] = img.getRadialDistortionCenterShift().y;
00186                 m_center[0] = img.getSize().width() / 2.0;
00187                 m_center[1] = img.getSize().height() / 2.0;
00188             }
00189         }
00190     };
00191     void SaveToImgs()
00192     {
00193         for (unsigned int i = 0; i < 3; i++)
00194         {
00195             HuginBase::SrcPanoImage img = m_pano.getSrcImage(i);
00196             img.setHFOV(m_hfov[i]);
00197             std::vector<double> radialDist(4);
00198             radialDist[0] = m_dist[i][0];
00199             radialDist[1] = m_dist[i][1];
00200             radialDist[2] = m_dist[i][2];
00201             radialDist[3] = 1 - radialDist[0] - radialDist[1] - radialDist[2];
00202             img.setRadialDistortion(radialDist);
00203             img.setRadialDistortionCenterShift(hugin_utils::FDiff2D(m_shift[0], m_shift[1]));
00204             m_pano.setSrcImage(i, img);
00205         }
00206     };
00207 };
00208 
00209 void get_optvars(HuginBase::OptimizeVector& _retval)
00210 {
00211     HuginBase::OptimizeVector optvars;
00212     std::set<std::string> vars = g_param.optvars;
00213     optvars.push_back(vars);
00214     optvars.push_back(std::set<std::string>());
00215     /* NOTE: delete "d" and "e" if they should be optimized,
00216     they are linked and always will be */
00217     vars.erase("d");
00218     vars.erase("e");
00219     optvars.push_back(vars);
00220     _retval = optvars;
00221 }
00222 
00223 // dummy panotools progress functions
00224 static int ptProgress(int command, char* argument)
00225 {
00226     return 1;
00227 }
00228 static int ptinfoDlg(int command, char* argument)
00229 {
00230     return 1;
00231 }
00232 
00233 // Method 0: using PTOptimizer
00234 //   PTOptimizer minimizes the tangential and sagittal distance between the points
00235 int optimize_old(HuginBase::Panorama& pano)
00236 {
00237     if (g_param.verbose == 0)
00238     {
00239         // deactive PTOptimizer status information if -v is not given
00240         PT_setProgressFcn(ptProgress);
00241         PT_setInfoDlgFcn(ptinfoDlg);
00242     };
00243 
00244     HuginBase::OptimizeVector optvars;
00245     get_optvars(optvars);
00246     pano.setOptimizeVector(optvars);
00247     HuginBase::PTools::optimize(pano);
00248     return 0;
00249 }
00250 
00251 inline double weightHuber(double x, double sigma)
00252 {
00253     if (fabs(x) > sigma)
00254     {
00255         x = sqrt(sigma*(2.0*fabs(x) - sigma));
00256     }
00257     return x;
00258 }
00259 
00260 void optGetError(double* p, double* x, int m, int n, void* data)
00261 {
00262     OptimData* dat = static_cast<OptimData*>(data);
00263     dat->FromX(p);
00264 
00265     /* compute new a,b,c,d from a,b,c,v */
00266     double dist[3][4];
00267     for (unsigned int i = 0; i<3; i++)
00268     {
00269         double scale = dat->m_hfov[1] / dat->m_hfov[i];
00270         for (unsigned int j = 0; j<3; j++)
00271         {
00272             dist[i][j] = dat->m_dist[i][j] * pow(scale, (int)(4 - j));
00273         }
00274         dist[i][3] = scale*(1 - dat->m_dist[i][0] - dat->m_dist[i][1] - dat->m_dist[i][2]);
00275     }
00276 
00277     double center[2];
00278     center[0] = dat->m_center[0] + dat->m_shift[0];
00279     center[1] = dat->m_center[1] + dat->m_shift[1];
00280 
00281     double base_size = std::min(dat->m_center[0], dat->m_center[1]);
00282 
00283     HuginBase::CPVector newCPs;
00284 
00285     unsigned int noPts = dat->m_pano.getNrOfCtrlPoints();
00286     // loop over all points to calculate the error
00287     for (unsigned int ptIdx = 0; ptIdx < noPts; ptIdx++)
00288     {
00289         const HuginBase::ControlPoint& cp = dat->m_pano.getCtrlPoint(ptIdx);
00290 
00291         double dist_p1 = vigra::hypot(cp.x1 - center[0], cp.y1 - center[1]);
00292         double dist_p2 = vigra::hypot(cp.x2 - center[0], cp.y2 - center[1]);
00293 
00294         if (cp.image1Nr == 1)
00295         {
00296             double base_dist = dist_p1 / base_size;
00297             double corr_dist_p1 = dist[cp.image2Nr][0] * pow(base_dist, 4) +
00298                 dist[cp.image2Nr][1] * pow(base_dist, 3) +
00299                 dist[cp.image2Nr][2] * pow(base_dist, 2) +
00300                 dist[cp.image2Nr][3] * base_dist;
00301             corr_dist_p1 *= base_size;
00302             x[ptIdx] = corr_dist_p1 - dist_p2;
00303         }
00304         else
00305         {
00306             double base_dist = dist_p2 / base_size;
00307             double corr_dist_p2 = dist[cp.image1Nr][0] * pow(base_dist, 4) +
00308                 dist[cp.image1Nr][1] * pow(base_dist, 3) +
00309                 dist[cp.image1Nr][2] * pow(base_dist, 2) +
00310                 dist[cp.image1Nr][3] * base_dist;
00311             corr_dist_p2 *= base_size;
00312             x[ptIdx] = corr_dist_p2 - dist_p1;
00313         }
00314 
00315         HuginBase::ControlPoint newcp = cp;
00316         newcp.error = fabs(x[ptIdx]);
00317         newCPs.push_back(newcp);
00318 
00319         dat->m_pano.getCtrlPoint(ptIdx);
00320         // use huber robust estimator
00321         if (dat->huberSigma > 0)
00322         {
00323             x[ptIdx] = weightHuber(x[ptIdx], dat->huberSigma);
00324         }
00325     }
00326 
00327     dat->m_pano.updateCtrlPointErrors(newCPs);
00328 }
00329 
00330 int optVis(double* p, double* x, int m, int n, int iter, double sqerror, void* data)
00331 {
00332     return 1;
00333     /*    OptimData * dat = (OptimData *) data;
00334     char tmp[200];
00335     tmp[199] = 0;
00336     double error = sqrt(sqerror/n)*255;
00337     snprintf(tmp,199, "Iteration: %d, error: %f", iter, error);
00338     return dat->m_progress.increaseProgress(0.0, tmp) ? 1 : 0 ; */
00339 }
00340 
00341 // Method 1: minimize only the center distance difference (sagittal distance) of the points
00342 //   the tangential distance is not of interest for TCA correction, 
00343 //   and is caused by the limited accuracy of the fine tune function, especially close the the edge of the fisheye image
00344 void optimize_new(HuginBase::PanoramaData& pano)
00345 {
00346     HuginBase::OptimizeVector optvars;
00347     get_optvars(optvars);
00348 
00349     int nMaxIter = 1000;
00350     OptimData data(pano, optvars, 0.5, nMaxIter);
00351 
00352     int ret;
00353     double info[LM_INFO_SZ];
00354 
00355     // parameters
00356     int m = data.m_mapping.size();
00357     vigra::ArrayVector<double> p(m, 0.0);
00358 
00359     // vector for errors
00360     int n = pano.getNrOfCtrlPoints();
00361     vigra::ArrayVector<double> x(n, 0.0);
00362 
00363     data.ToX(p.begin());
00364     if (g_param.verbose > 0)
00365     {
00366         fprintf(stderr, "Parameters before optimization: ");
00367         for (int i = 0; i<m; ++i)
00368         {
00369             fprintf(stderr, "%.7g ", p[i]);
00370         }
00371         fprintf(stderr, "\n");
00372     }
00373 
00374     // TODO: setup optimization options with some good defaults.
00375     double optimOpts[5];
00376 
00377     optimOpts[0] = 1e-5;   // init mu
00378     // stop thresholds
00379     optimOpts[1] = 1e-7;   // ||J^T e||_inf
00380     optimOpts[2] = 1e-10;   // ||Dp||_2
00381     optimOpts[3] = 1e-3;   // ||e||_2
00382     // difference mode
00383     optimOpts[4] = LM_DIFF_DELTA;
00384 
00385     //    data.huberSigma = 0;
00386 
00387     ret = dlevmar_dif(&optGetError, &optVis, &(p[0]), &(x[0]), m, n, nMaxIter, NULL, info, NULL, NULL, &data);  // no jacobian
00388     // copy to source images (data.m_imgs)
00389     data.SaveToImgs();
00390     // calculate error at solution
00391     data.huberSigma = 0;
00392     optGetError(&(p[0]), &(x[0]), m, n, &data);
00393     double error = 0;
00394     for (int i = 0; i<n; i++)
00395     {
00396         error += x[i] * x[i];
00397     }
00398     error = sqrt(error / n);
00399 
00400     if (g_param.verbose > 0)
00401     {
00402         fprintf(stderr, "Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
00403         for (int i = 0; i<m; ++i)
00404         {
00405             fprintf(stderr, "%.7g ", p[i]);
00406         }
00407         fprintf(stderr, "\n\nMinimization info:\n");
00408         for (int i = 0; i<LM_INFO_SZ; ++i)
00409         {
00410             fprintf(stderr, "%g ", info[i]);
00411         }
00412         fprintf(stderr, "\n");
00413     }
00414 }
00415 
00416 static void usage(const char* name)
00417 {
00418     std::cerr << name << ": Parameter estimation of transverse chromatic abberations" << std::endl
00419          << name << " version " << hugin_utils::GetHuginVersion() << std::endl
00420          << std::endl
00421          << "Usage: " << name  << " [options] <inputfile>" << std::endl
00422          << "  option are: " << std::endl
00423          << "    -h            Display help (this text)" << std::endl
00424          << "    -l            input file is PTO file instead of image" << std::endl
00425          << "    -m method     optimization method (0 normal, 1 newfit)" << std::endl
00426          << "    -o optvars    string of variables to optimize (\"abcvde\")" << std::endl
00427          << "    -r            reset values (this will zero a,b,c,d,e params and set v to 10)" << std::endl
00428          << "                  makes sense only with -l option" << std::endl
00429          << "    -s <scale>    Scale for corner detection" << std::endl
00430          << "    -n <number>   number of points per grid cell (default: 10)" << std::endl
00431          << "    -g <number>   divide image in <number>x<number> grid cells (default: 10)" << std::endl
00432          << "    -t num        Remove all control points with an error higher than num pixels (default: 1.5)" << std::endl
00433          << "    -v            Verbose" << std::endl
00434          << "    --save-into-database  Saves the tca data into Hugin lens database" << std::endl
00435          << "    -w filename   write PTO file" << std::endl
00436          << "    -R <r>        Use this file as red channel" << std::endl
00437          << "    -G <g>        Use this file as green channel" << std::endl
00438          << "    -B <b>        Use this file as blue channel" << std::endl
00439          << std::endl
00440          << "  <inputfile> is the base name of 4 image files:" << std::endl
00441          << "    <inputfile>        Colour file to compute TCA parameters" << std::endl
00442          << "    red_<inputfile>    Red channel of <inputfile>" << std::endl
00443          << "    green_<inputfile>  Green channel of <inputfile>" << std::endl
00444          << "    blue_<inputfile>   Blue channel of <inputfile>" << std::endl
00445          << "    The channel images must be colour images with 3 identical channels." << std::endl
00446          << "    If any of -R, -G, or -B is given, this file name is used instead of the derived name." << std::endl
00447          << std::endl
00448          << "  Output:" << std::endl
00449          << "    commandline arguments for fulla" << std::endl;
00450 }
00451 
00452 static hugin_omp::Lock lock;
00453 typedef std::multimap<double, vigra::Diff2D> MapPoints;
00454 
00455 template <class ImageType>
00456 void createCtrlPoints(HuginBase::Panorama& pano, const ImageType& img, int imgRedNr, int imgGreenNr, int imgBlueNr, double scale, int nPoints, unsigned grid)
00457 {
00458     vigra::BasicImage<vigra::RGBValue<vigra::UInt8> > img8(img.size());
00459 
00460     double ratio = 255.0/vigra_ext::LUTTraits<typename ImageType::value_type>::max();
00461     transformImage(srcImageRange(img), destImage(img8),
00462         vigra::functor::Arg1()*vigra::functor::Param(ratio));
00463     if (g_param.verbose > 0)
00464     {
00465         std::cout << "image8 size:" << img8.size() << std::endl;
00466     };
00468     // find interesting corners using harris corner detector
00469 
00470     if (g_param.verbose > 0)
00471     {
00472         std::cout << "Finding control points... " << std::endl;
00473     }
00474     const long templWidth = 29;
00475     const long sWidth = 29 + 11;
00476 
00477 
00478     vigra::Size2D size(img8.width(), img8.height());
00479     std::vector<vigra::Rect2D> rects;
00480     for (unsigned party = 0; party < grid; party++)
00481     {
00482         for (unsigned partx = 0; partx < grid; partx++)
00483         {
00484             // run corner detector only in current sub-region (saves a lot of memory for big images)
00485             vigra::Rect2D rect(partx*size.x / grid, party*size.y / grid,
00486                 (partx + 1)*size.x / grid, (party + 1)*size.y / grid);
00487             rect &= vigra::Rect2D(size);
00488             if (rect.width()>0 && rect.height()>0)
00489             {
00490                 rects.push_back(rect);
00491             };
00492         };
00493     };
00494 
00495     #pragma omp parallel for schedule(dynamic)
00496     for (int i = 0; i < rects.size(); ++i)
00497     {
00498         MapPoints points;
00499         vigra::Rect2D rect(rects[i]);
00500         vigra_ext::findInterestPointsPartial(vigra::srcImageRange(img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >()), rect, scale, 5 * nPoints, points);
00501 
00502         // loop over all points, starting with the highest corner score
00503         HuginBase::CPVector cps;
00504         size_t nBad = 0;
00505         for (MapPoints::const_reverse_iterator it = points.rbegin(); it != points.rend(); ++it)
00506         {
00507             if (cps.size() >= nPoints)
00508             {
00509                 // we have enough points, stop
00510                 break;
00511             };
00512 
00513             // Green <-> Red
00514             HuginBase::ControlPoint p1(imgGreenNr, it->second.x, it->second.y, imgRedNr, it->second.x, it->second.y);
00515             vigra::Diff2D roundP1(hugin_utils::roundi(p1.x1), hugin_utils::roundi(p1.y1));
00516             vigra::Diff2D roundP2(hugin_utils::roundi(p1.x2), hugin_utils::roundi(p1.y2));
00517 
00518             vigra_ext::CorrelationResult res = vigra_ext::PointFineTune(
00519                 img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >(),
00520                 roundP1, templWidth,
00521                 img8, vigra::RedAccessor<vigra::RGBValue<vigra::UInt8> >(),
00522                 roundP2, sWidth);
00523 
00524             if (res.maxi > 0.98)
00525             {
00526                 p1.x1 = roundP1.x;
00527                 p1.y1 = roundP1.y;
00528                 p1.x2 = res.maxpos.x;
00529                 p1.y2 = res.maxpos.y;
00530                 p1.error = res.maxi;
00531                 cps.push_back(p1);
00532             }
00533             else
00534             {
00535                 ++nBad;
00536             };
00537 
00538             // Green <-> Blue
00539             HuginBase::ControlPoint p2(imgGreenNr, it->second.x, it->second.y, imgBlueNr, it->second.x, it->second.y);
00540             roundP1 = vigra::Diff2D(hugin_utils::roundi(p2.x1), hugin_utils::roundi(p2.y1));
00541             roundP2 = vigra::Diff2D(hugin_utils::roundi(p2.x2), hugin_utils::roundi(p2.y2));
00542 
00543             res = vigra_ext::PointFineTune(
00544                 img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >(), roundP1, templWidth,
00545                 img8, vigra::BlueAccessor<vigra::RGBValue<vigra::UInt8> >(), roundP2, sWidth);
00546 
00547             if (res.maxi > 0.98)
00548             {
00549                 p2.x1 = roundP1.x;
00550                 p2.y1 = roundP1.y;
00551                 p2.x2 = res.maxpos.x;
00552                 p2.y2 = res.maxpos.y;
00553                 p2.error = res.maxi;
00554                 cps.push_back(p2);
00555             }
00556             else
00557             {
00558                 ++nBad;
00559             };
00560         }
00561         if (g_param.verbose > 0)
00562         {
00563             std::ostringstream buf;
00564             buf << "Number of good matches: " << cps.size() << ", bad matches: " << nBad << std::endl;
00565             std::cout << buf.str();
00566         }
00567         if (!cps.empty())
00568         {
00569             hugin_omp::ScopedLock sl(lock);
00570             for (HuginBase::CPVector::const_iterator it = cps.begin(); it != cps.end(); ++it)
00571             {
00572                 pano.addCtrlPoint(*it);
00573             };
00574         };
00575     };
00576 };
00577 
00578 int main2(HuginBase::Panorama& pano);
00579 
00580 template <class PixelType>
00581 int processImg(const char* filename)
00582 {
00583     typedef vigra::BasicImage<PixelType> ImageType;
00584     try
00585     {
00586         // load first image
00587         vigra::ImageImportInfo imgInfo(filename);
00588         ImageType imgOrig(imgInfo.size());
00589         
00590         const int bands = imgInfo.numBands();
00591         const int extraBands = imgInfo.numExtraBands();
00592 
00593         if (!(bands == 3 || (bands == 4 && extraBands == 1)))
00594         {
00595             std::cerr << "Unsupported number of bands!";
00596             exit(-1);
00597         }
00598 
00599         if (bands == 3)
00600         {
00601             vigra::importImage(imgInfo, vigra::destImage(imgOrig));
00602         }
00603         else
00604         {
00605             vigra::BImage alpha(imgInfo.size());
00606             vigra::importImageAlpha(imgInfo, destImage(imgOrig), destImage(alpha));
00607         };
00608         HuginBase::Panorama pano;
00609         // add the first image to the panorama object
00610         HuginBase::StandardImageVariableGroups variable_groups(pano);
00611         HuginBase::ImageVariableGroup& lenses = variable_groups.getLenses();
00612 
00613         std::string red_name;
00614         if( g_param.red_name.size())
00615         {
00616             red_name=g_param.red_name;
00617         }
00618         else
00619         {
00620             red_name=std::string("red_")+filename;
00621         }
00622 
00623         HuginBase::SrcPanoImage srcRedImg;
00624         srcRedImg.setSize(imgInfo.size());
00625         srcRedImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00626         srcRedImg.setHFOV(10);
00627         srcRedImg.setCropFactor(1);
00628         srcRedImg.setFilename(red_name);
00629         int imgRedNr = pano.addImage(srcRedImg);
00630         lenses.updatePartNumbers();
00631         lenses.switchParts(imgRedNr, 0);
00632 
00633         std::string green_name;
00634         if( g_param.green_name.size())
00635         {
00636             green_name=g_param.green_name;
00637         }
00638         else
00639         {
00640             green_name=std::string("green_")+filename;
00641         }
00642 
00643         HuginBase::SrcPanoImage srcGreenImg;
00644         srcGreenImg.setSize(imgInfo.size());
00645         srcGreenImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00646         srcGreenImg.setHFOV(10);
00647         srcGreenImg.setCropFactor(1);
00648         srcGreenImg.setFilename(green_name);
00649         int imgGreenNr = pano.addImage(srcGreenImg);
00650         lenses.updatePartNumbers();
00651         lenses.switchParts(imgGreenNr, 0);
00652 
00653         std::string blue_name;
00654         if( g_param.blue_name.size())
00655         {
00656             blue_name=g_param.blue_name;
00657         }
00658         else
00659         {
00660             blue_name=std::string("blue_")+filename;
00661         }
00662 
00663         HuginBase::SrcPanoImage srcBlueImg;
00664         srcBlueImg.setSize(imgInfo.size());
00665         srcBlueImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00666         srcBlueImg.setHFOV(10);
00667         srcBlueImg.setCropFactor(1);
00668         srcBlueImg.setFilename(blue_name);
00669         int imgBlueNr = pano.addImage(srcBlueImg);
00670         lenses.updatePartNumbers();
00671         lenses.switchParts(imgBlueNr, 0);
00672 
00673         // lens variables are linked by default. Unlink the field of view and
00674         // the radial distortion.
00675         lenses.unlinkVariablePart(HuginBase::ImageVariableGroup::IVE_HFOV, 0);
00676         lenses.unlinkVariablePart(HuginBase::ImageVariableGroup::IVE_RadialDistortion, 0);
00677 
00678         // setup output to be exactly similar to input image
00679         HuginBase::PanoramaOptions opts;
00680         opts.setProjection(HuginBase::PanoramaOptions::RECTILINEAR);
00681         opts.setHFOV(srcGreenImg.getHFOV(), false);
00682         opts.setWidth(srcGreenImg.getSize().x, false);
00683         opts.setHeight(srcGreenImg.getSize().y);
00684 
00685         // output to tiff format
00686         opts.outputFormat = HuginBase::PanoramaOptions::TIFF_m;
00687         opts.tiff_saveROI = false;
00688         // m estimator, to be more robust against points on moving objects
00689         opts.huberSigma = 0.5;
00690         pano.setOptions(opts);
00691 
00692         createCtrlPoints(pano, imgOrig, imgRedNr, imgGreenNr, imgBlueNr, g_param.scale, g_param.nPoints, g_param.grid);
00693 
00694         main2(pano);
00695     }
00696     catch (std::exception& e)
00697     {
00698         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
00699         return 1;
00700     }
00701     return 0;
00702 }
00703 
00704 int processPTO(const char* filename)
00705 {
00706     HuginBase::Panorama pano;
00707 
00708     std::ifstream ptofile(filename);
00709     if (ptofile.bad())
00710     {
00711         std::cerr << "could not open script : " << filename << std::endl;
00712         return 1;
00713     }
00714     pano.setFilePrefix(hugin_utils::getPathPrefix(filename));
00715     AppBase::DocumentData::ReadWriteError err = pano.readData(ptofile);
00716     if (err != AppBase::DocumentData::SUCCESSFUL)
00717     {
00718         std::cerr << "error while parsing script: " << filename << std::endl;
00719         return 1;
00720     }
00721 
00722     return main2(pano);
00723 }
00724 
00725 void resetValues(HuginBase::Panorama& pano)
00726 {
00727     for (unsigned int i=0; i < 3; i++)
00728     {
00729         HuginBase::SrcPanoImage img = pano.getSrcImage(i);
00730 
00731         img.setHFOV(10);
00732 
00733         std::vector<double> radialDist(4);
00734         radialDist[0] = 0;
00735         radialDist[1] = 0;
00736         radialDist[2] = 0;
00737         radialDist[3] = 1;
00738         img.setRadialDistortion(radialDist);
00739 
00740         img.setRadialDistortionCenterShift(hugin_utils::FDiff2D(0, 0));
00741 
00742         pano.setSrcImage(i, img);
00743     }
00744 }
00745 
00746 void print_result(HuginBase::Panorama& pano)
00747 {
00748     double dist[3][3]; // a,b,c for all imgs
00749     double shift[2];   // x,y shift
00750     double hfov[3];
00751 
00752     for (unsigned int i=0; i < 3; i++)
00753     {
00754         HuginBase::SrcPanoImage img = pano.getSrcImage(i);
00755         hfov[i] = img.getHFOV();
00756         dist[i][0] = img.getRadialDistortion()[0];
00757         dist[i][1] = img.getRadialDistortion()[1];
00758         dist[i][2] = img.getRadialDistortion()[2];
00759         if (i == 0)
00760         {
00761             shift[0] = img.getRadialDistortionCenterShift().x;
00762             shift[1] = img.getRadialDistortionCenterShift().y;
00763         }
00764     }
00765 
00766     /* compute new a,b,c,d from a,b,c,v */
00767     double distnew[3][4];
00768     for (unsigned int i=0 ; i<3 ; i++)
00769     {
00770         double scale = hfov[1] / hfov[i];
00771         for (unsigned int j=0 ; j<3 ; j++)
00772         {
00773             distnew[i][j] = dist[i][j]*pow(scale, (int)(4-j));
00774         }
00775         distnew[i][3] = scale*(1 - dist[i][0] - dist[i][1] - dist[i][2]);
00776     }
00777 
00778     if (hugin_utils::roundi(shift[0]) == 0 && hugin_utils::roundi(shift[1]) == 0)
00779     {
00780         fprintf(stdout, "-r %.7f:%.7f:%.7f:%.7f "
00781             "-b %.7f:%.7f:%.7f:%.7f ",
00782             distnew[0][0], distnew[0][1], distnew[0][2], distnew[0][3],
00783             distnew[2][0], distnew[2][1], distnew[2][2], distnew[2][3]);
00784     }
00785     else
00786     {
00787         fprintf(stdout, "-r %.7f:%.7f:%.7f:%.7f "
00788             "-b %.7f:%.7f:%.7f:%.7f "
00789             "-x %d:%d\n",
00790             distnew[0][0], distnew[0][1], distnew[0][2], distnew[0][3],
00791             distnew[2][0], distnew[2][1], distnew[2][2], distnew[2][3],
00792             hugin_utils::roundi(shift[0]), hugin_utils::roundi(shift[1]));
00793     };
00794     if (g_param.saveDB)
00795     {
00796         // save information into database
00797         // read EXIF data, using routines of HuginBase::SrcPanoImage
00798         HuginBase::SrcPanoImage img;
00799         img.setFilename(g_param.basename);
00800         img.readEXIF();
00801         std::string lensname = img.getDBLensName();
00802         if (lensname.empty())
00803         {
00804             // no suitable lensname found in exif data, ask user
00805             std::cout << std::endl << "For saving information tca data into database no suitable information" << std::endl
00806                 << "found in EXIF data." << std::endl
00807                 << "For fixed lens cameras leave lensname empty." << std::endl << std::endl << "Lensname: ";
00808             char input[256];
00809             std::cin.getline(input, 255);
00810             lensname = hugin_utils::StrTrim(std::string(input));
00811             if (lensname.empty())
00812             {
00813                 std::string camMaker;
00814                 std::string camModel;
00815                 std::cout << "Camera maker: ";
00816                 while (camMaker.empty())
00817                 {
00818                     std::cin.getline(input, 255);
00819                     camMaker = hugin_utils::StrTrim(std::string(input));
00820                 };
00821                 std::cout << "Camera model: ";
00822                 while (camModel.empty())
00823                 {
00824                     std::cin.getline(input, 255);
00825                     camModel = hugin_utils::StrTrim(std::string(input));
00826                 };
00827                 lensname = camMaker.append("|").append(camModel);
00828             };
00829         };
00830         double focal = img.getExifFocalLength();
00831         if (fabs(focal) < 0.1f)
00832         {
00833             std::cout << "Real focal length (in mm): ";
00834             while (fabs(focal) < 0.1f)
00835             {
00836                 std::cin >> focal;
00837                 focal = fabs(focal);
00838             };
00839         };
00840         HuginBase::LensDB::LensDB& lensDB = HuginBase::LensDB::LensDB::GetSingleton();
00841         std::vector<double> redDistData(4, 0.0f);
00842         std::vector<double> blueDistData(4, 0.0f);
00843         redDistData[0] = distnew[0][0];
00844         redDistData[1] = distnew[0][1];
00845         redDistData[2] = distnew[0][2];
00846         redDistData[3] = distnew[0][3];
00847         blueDistData[0] = distnew[2][0];
00848         blueDistData[1] = distnew[2][1];
00849         blueDistData[2] = distnew[2][2];
00850         blueDistData[3] = distnew[2][3];
00851         if (lensDB.SaveTCA(lensname, focal, redDistData, blueDistData))
00852         {
00853             std::cout << std::endl << std::endl << "TCA data for " << lensname << " @ " << focal << " mm successful saved into database." << std::endl;
00854         }
00855         else
00856         {
00857             std::cout << std::endl << std::endl << "Could not save data into database." << std::endl;
00858         };
00859         HuginBase::LensDB::LensDB::Clean();
00860     }
00861 }
00862 
00863 int main2(HuginBase::Panorama& pano)
00864 {
00865     if (g_param.reset)
00866     {
00867         resetValues(pano);
00868     }
00869 
00870     for (int i = 0; i < 10; i++)
00871     {
00872         if (g_param.optMethod == 0)
00873         {
00874             optimize_old(pano);
00875         }
00876         else if (g_param.optMethod == 1)
00877         {
00878             optimize_new(pano);
00879         }
00880         else
00881         {
00882             std::cerr << "Unknown optimizer strategy." << std::endl
00883                 << "Using newfit method." << std::endl;
00884             g_param.optMethod = 1;
00885             optimize_new(pano);
00886         };
00887 
00888         HuginBase::CPVector cps = pano.getCtrlPoints();
00889         HuginBase::CPVector newCPs;
00890         for (HuginBase::CPVector::const_iterator it = cps.begin(); it != cps.end(); ++it)
00891         {
00892             if (it->error < g_param.cpErrorThreshold)
00893             {
00894                 newCPs.push_back(*it);
00895             }
00896         }
00897         if (g_param.verbose > 0)
00898         {
00899             std::cerr << "Ctrl points before pruning: " << cps.size() << ", after: " << newCPs.size() << std::endl;
00900         }
00901         pano.setCtrlPoints(newCPs);
00902 
00903         if (cps.size() == newCPs.size())
00904         {
00905             // no points were removed, do not re-optimize
00906             break;
00907         }
00908     }
00909 
00910     if (!g_param.ptoOutputFile.empty())
00911     {
00912         HuginBase::OptimizeVector optvars;
00913         get_optvars(optvars);
00914         HuginBase::UIntSet allImgs;
00915         fill_set(allImgs, 0, pano.getNrOfImages() - 1);
00916         std::ofstream script(g_param.ptoOutputFile.c_str());
00917         pano.printPanoramaScript(script, optvars, pano.getOptions(), allImgs, true, "");
00918     }
00919 
00920     print_result(pano);
00921     return 0;
00922 }
00923 
00924 int main(int argc, char* argv[])
00925 {
00926     // parse arguments
00927     const char* optstring = "hlm:o:rt:vw:R:G:B:s:g:n:";
00928     enum
00929     {
00930         SAVEDATABASE=1000
00931     };
00932     static struct option longOptions[] =
00933     {
00934         { "save-into-database", no_argument, NULL, SAVEDATABASE },
00935         { "help", no_argument, NULL, 'h' },
00936         0
00937     };
00938     int c;
00939     bool parameter_request_seen=false;
00940     int optionIndex = 0;
00941     opterr = 0;
00942 
00943     while ((c = getopt_long(argc, argv, optstring, longOptions, &optionIndex)) != -1)
00944     {
00945         switch (c)
00946         {
00947             case 'h':
00948                 usage(hugin_utils::stripPath(argv[0]).c_str());
00949                 return 0;
00950             case 'l':
00951                 g_param.load = true;
00952                 break;
00953             case 'm':
00954                 g_param.optMethod = atoi(optarg);
00955                 break;
00956             case 'o':
00957                 {
00958                     char* optptr = optarg;
00959                     while (*optptr != 0)
00960                     {
00961                         if ((*optptr == 'a') || (*optptr == 'b') ||
00962                                 (*optptr == 'c') || (*optptr == 'v') ||
00963                                 (*optptr == 'd') || (*optptr == 'e'))
00964                         {
00965                             g_param.optvars.insert(std::string(optptr, 1));
00966                         }
00967                         optptr++;
00968                     }
00969                     parameter_request_seen=true;
00970                 }
00971                 break;
00972             case 'r':
00973                 g_param.reset = true;
00974                 break;
00975             case 't':
00976                 g_param.cpErrorThreshold = atof(optarg);
00977                 if (g_param.cpErrorThreshold <= 0)
00978                 {
00979                     std::cerr << "Invalid parameter: control point error threshold (-t) must be greater than 0" << std::endl;
00980                     return 1;
00981                 }
00982                 break;
00983             case 'v':
00984                 g_param.verbose++;
00985                 break;
00986             case 'w':
00987                 g_param.ptoOutputFile = optarg;
00988                 break;
00989             case 'R':
00990                 g_param.red_name = optarg;
00991                 break;
00992             case 'G':
00993                 g_param.green_name = optarg;
00994                 break;
00995             case 'B':
00996                 g_param.blue_name = optarg;
00997                 break;
00998             case 's':
00999                 g_param.scale=atof( optarg);
01000                 break;
01001             case 'n':
01002                 g_param.nPoints=atoi( optarg);
01003                 break;
01004             case 'g':
01005                 g_param.grid=atoi(optarg);
01006                 break;
01007             case SAVEDATABASE:
01008                 g_param.saveDB = true;
01009                 break;
01010             default:
01011                 std::cerr << "Invalid parameter: '" << argv[optind-1] << " " << optarg << "'" << std::endl;
01012                 usage(hugin_utils::stripPath(argv[0]).c_str());
01013                 return 1;
01014         }
01015     };
01016 
01017     if ((argc - optind) != 1)
01018     {
01019         usage(hugin_utils::stripPath(argv[0]).c_str());
01020         return 1;
01021     }
01022 
01023     // If no parameters were requested to be optimised, we optimize the
01024     // default parameters.
01025     if ( !parameter_request_seen)
01026     {
01027         for ( const char* dop=DEFAULT_OPTIMISATION_PARAMETER;
01028                 *dop != 0; ++dop)
01029         {
01030             g_param.optvars.insert( std::string( dop, 1));
01031         }
01032     }
01033 
01034     // Program will crash if nothing is to be optimised.
01035     if ( g_param.optvars.empty())
01036     {
01037         std::cerr << "No parameters to optimize." << std::endl;
01038         usage(hugin_utils::stripPath(argv[0]).c_str());
01039         return 1;
01040     }
01041 
01042     if (!g_param.load)
01043     {
01044         g_param.basename = argv[optind];
01045         vigra::ImageImportInfo firstImgInfo(argv[optind]);
01046         std::string pixelType = firstImgInfo.getPixelType();
01047         if (pixelType == "UINT8")
01048         {
01049             return processImg<vigra::RGBValue<vigra::UInt8> >(argv[optind]);
01050         }
01051         else if (pixelType == "INT16")
01052         {
01053             return processImg<vigra::RGBValue<vigra::Int16> >(argv[optind]);
01054         }
01055         else if (pixelType == "UINT16")
01056         {
01057             return processImg<vigra::RGBValue<vigra::UInt16> >(argv[optind]);
01058         }
01059         else
01060         {
01061             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01062             return 1;
01063         }
01064     }
01065     else
01066     {
01067         return processPTO(argv[optind]);
01068     }
01069 
01070     return 0;
01071 }
01072 

Generated on 7 Feb 2016 for Hugintrunk by  doxygen 1.4.7