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

Generated on 4 Sep 2015 for Hugintrunk by  doxygen 1.4.7