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     ret = dlevmar_dif(&optGetError, &optVis, &(p[0]), &(x[0]), m, n, nMaxIter, NULL, info, NULL, NULL, &data);  // no jacobian
00375     // copy to source images (data.m_imgs)
00376     data.SaveToImgs();
00377     // calculate error at solution
00378     data.huberSigma = 0;
00379     optGetError(&(p[0]), &(x[0]), m, n, &data);
00380     double error = 0;
00381     for (int i = 0; i<n; i++)
00382     {
00383         error += x[i] * x[i];
00384     }
00385     error = sqrt(error / n);
00386 
00387     if (g_param.verbose > 0)
00388     {
00389         fprintf(stderr, "Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
00390         for (int i = 0; i<m; ++i)
00391         {
00392             fprintf(stderr, "%.7g ", p[i]);
00393         }
00394         fprintf(stderr, "\n\nMinimization info:\n");
00395         for (int i = 0; i<LM_INFO_SZ; ++i)
00396         {
00397             fprintf(stderr, "%g ", info[i]);
00398         }
00399         fprintf(stderr, "\n");
00400     }
00401 }
00402 
00403 static void usage(const char* name)
00404 {
00405     std::cerr << name << ": Parameter estimation of transverse chromatic abberations" << std::endl
00406          << name << " version " << hugin_utils::GetHuginVersion() << std::endl
00407          << std::endl
00408          << "Usage: " << name  << " [options] <inputfile>" << std::endl
00409          << "  option are: " << std::endl
00410          << "    -h            Display help (this text)" << std::endl
00411          << "    -l            input file is PTO file instead of image" << std::endl
00412          << "    -m method     optimization method (0 normal, 1 newfit)" << std::endl
00413          << "    -o optvars    string of variables to optimize (\"abcvde\")" << std::endl
00414          << "    -r            reset values (this will zero a,b,c,d,e params and set v to 10)" << std::endl
00415          << "                  makes sense only with -l option" << std::endl
00416          << "    -s <scale>    Scale for corner detection" << std::endl
00417          << "    -n <number>   number of points per grid cell (default: 10)" << std::endl
00418          << "    -g <number>   divide image in <number>x<number> grid cells (default: 10)" << std::endl
00419          << "    -t num        Remove all control points with an error higher than num pixels (default: 1.5)" << std::endl
00420          << "    -v            Verbose" << std::endl
00421          << "    --save-into-database  Saves the tca data into Hugin lens database" << std::endl
00422          << "    -w filename   write PTO file" << std::endl
00423          << "    -R <r>        Use this file as red channel" << std::endl
00424          << "    -G <g>        Use this file as green channel" << std::endl
00425          << "    -B <b>        Use this file as blue channel" << std::endl
00426          << std::endl
00427          << "  <inputfile> is the base name of 4 image files:" << std::endl
00428          << "    <inputfile>        Colour file to compute TCA parameters" << std::endl
00429          << "    red_<inputfile>    Red channel of <inputfile>" << std::endl
00430          << "    green_<inputfile>  Green channel of <inputfile>" << std::endl
00431          << "    blue_<inputfile>   Blue channel of <inputfile>" << std::endl
00432          << "    The channel images must be colour images with 3 identical channels." << std::endl
00433          << "    If any of -R, -G, or -B is given, this file name is used instead of the derived name." << std::endl
00434          << std::endl
00435          << "  Output:" << std::endl
00436          << "    commandline arguments for fulla" << std::endl;
00437 }
00438 
00439 static hugin_omp::Lock lock;
00440 typedef std::multimap<double, vigra::Diff2D> MapPoints;
00441 
00442 template <class ImageType>
00443 void createCtrlPoints(HuginBase::Panorama& pano, const ImageType& img, int imgRedNr, int imgGreenNr, int imgBlueNr, double scale, int nPoints, unsigned grid)
00444 {
00445     vigra::BasicImage<vigra::RGBValue<vigra::UInt8> > img8(img.size());
00446 
00447     double ratio = 255.0/vigra_ext::LUTTraits<typename ImageType::value_type>::max();
00448     transformImage(srcImageRange(img), destImage(img8),
00449         vigra::functor::Arg1()*vigra::functor::Param(ratio));
00450     if (g_param.verbose > 0)
00451     {
00452         std::cout << "image8 size:" << img8.size() << std::endl;
00453     };
00455     // find interesting corners using harris corner detector
00456 
00457     if (g_param.verbose > 0)
00458     {
00459         std::cout << "Finding control points... " << std::endl;
00460     }
00461     const long templWidth = 29;
00462     const long sWidth = 29 + 11;
00463 
00464 
00465     vigra::Size2D size(img8.width(), img8.height());
00466     std::vector<vigra::Rect2D> rects;
00467     for (unsigned party = 0; party < grid; party++)
00468     {
00469         for (unsigned partx = 0; partx < grid; partx++)
00470         {
00471             // run corner detector only in current sub-region (saves a lot of memory for big images)
00472             vigra::Rect2D rect(partx*size.x / grid, party*size.y / grid,
00473                 (partx + 1)*size.x / grid, (party + 1)*size.y / grid);
00474             rect &= vigra::Rect2D(size);
00475             if (rect.width()>0 && rect.height()>0)
00476             {
00477                 rects.push_back(rect);
00478             };
00479         };
00480     };
00481 
00482     #pragma omp parallel for schedule(dynamic)
00483     for (int i = 0; i < rects.size(); ++i)
00484     {
00485         MapPoints points;
00486         vigra::Rect2D rect(rects[i]);
00487         vigra_ext::findInterestPointsPartial(vigra::srcImageRange(img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >()), rect, scale, 5 * nPoints, points);
00488 
00489         // loop over all points, starting with the highest corner score
00490         HuginBase::CPVector cps;
00491         size_t nBad = 0;
00492         for (MapPoints::const_reverse_iterator it = points.rbegin(); it != points.rend(); ++it)
00493         {
00494             if (cps.size() >= nPoints)
00495             {
00496                 // we have enough points, stop
00497                 break;
00498             };
00499 
00500             // Green <-> Red
00501             HuginBase::ControlPoint p1(imgGreenNr, it->second.x, it->second.y, imgRedNr, it->second.x, it->second.y);
00502             vigra::Diff2D roundP1(hugin_utils::roundi(p1.x1), hugin_utils::roundi(p1.y1));
00503             vigra::Diff2D roundP2(hugin_utils::roundi(p1.x2), hugin_utils::roundi(p1.y2));
00504 
00505             vigra_ext::CorrelationResult res = vigra_ext::PointFineTune(
00506                 img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >(),
00507                 roundP1, templWidth,
00508                 img8, vigra::RedAccessor<vigra::RGBValue<vigra::UInt8> >(),
00509                 roundP2, sWidth);
00510 
00511             if (res.maxi > 0.98)
00512             {
00513                 p1.x1 = roundP1.x;
00514                 p1.y1 = roundP1.y;
00515                 p1.x2 = res.maxpos.x;
00516                 p1.y2 = res.maxpos.y;
00517                 p1.error = res.maxi;
00518                 cps.push_back(p1);
00519             }
00520             else
00521             {
00522                 ++nBad;
00523             };
00524 
00525             // Green <-> Blue
00526             HuginBase::ControlPoint p2(imgGreenNr, it->second.x, it->second.y, imgBlueNr, it->second.x, it->second.y);
00527             roundP1 = vigra::Diff2D(hugin_utils::roundi(p2.x1), hugin_utils::roundi(p2.y1));
00528             roundP2 = vigra::Diff2D(hugin_utils::roundi(p2.x2), hugin_utils::roundi(p2.y2));
00529 
00530             res = vigra_ext::PointFineTune(
00531                 img8, vigra::GreenAccessor<vigra::RGBValue<vigra::UInt8> >(), roundP1, templWidth,
00532                 img8, vigra::BlueAccessor<vigra::RGBValue<vigra::UInt8> >(), roundP2, sWidth);
00533 
00534             if (res.maxi > 0.98)
00535             {
00536                 p2.x1 = roundP1.x;
00537                 p2.y1 = roundP1.y;
00538                 p2.x2 = res.maxpos.x;
00539                 p2.y2 = res.maxpos.y;
00540                 p2.error = res.maxi;
00541                 cps.push_back(p2);
00542             }
00543             else
00544             {
00545                 ++nBad;
00546             };
00547         }
00548         if (g_param.verbose > 0)
00549         {
00550             std::ostringstream buf;
00551             buf << "Number of good matches: " << cps.size() << ", bad matches: " << nBad << std::endl;
00552             std::cout << buf.str();
00553         }
00554         if (!cps.empty())
00555         {
00556             hugin_omp::ScopedLock sl(lock);
00557             for (HuginBase::CPVector::const_iterator it = cps.begin(); it != cps.end(); ++it)
00558             {
00559                 pano.addCtrlPoint(*it);
00560             };
00561         };
00562     };
00563 };
00564 
00565 int main2(HuginBase::Panorama& pano);
00566 
00567 template <class PixelType>
00568 int processImg(const char* filename)
00569 {
00570     typedef vigra::BasicImage<PixelType> ImageType;
00571     try
00572     {
00573         // load first image
00574         vigra::ImageImportInfo imgInfo(filename);
00575         ImageType imgOrig(imgInfo.size());
00576         
00577         const int bands = imgInfo.numBands();
00578         const int extraBands = imgInfo.numExtraBands();
00579 
00580         if (!(bands == 3 || (bands == 4 && extraBands == 1)))
00581         {
00582             std::cerr << "Unsupported number of bands!";
00583             exit(-1);
00584         }
00585 
00586         if (bands == 3)
00587         {
00588             vigra::importImage(imgInfo, vigra::destImage(imgOrig));
00589         }
00590         else
00591         {
00592             vigra::BImage alpha(imgInfo.size());
00593             vigra::importImageAlpha(imgInfo, destImage(imgOrig), destImage(alpha));
00594         };
00595         HuginBase::Panorama pano;
00596         // add the first image to the panorama object
00597         HuginBase::StandardImageVariableGroups variable_groups(pano);
00598         HuginBase::ImageVariableGroup& lenses = variable_groups.getLenses();
00599 
00600         std::string red_name;
00601         if( g_param.red_name.size())
00602         {
00603             red_name=g_param.red_name;
00604         }
00605         else
00606         {
00607             red_name=std::string("red_")+filename;
00608         }
00609 
00610         HuginBase::SrcPanoImage srcRedImg;
00611         srcRedImg.setSize(imgInfo.size());
00612         srcRedImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00613         srcRedImg.setHFOV(10);
00614         srcRedImg.setCropFactor(1);
00615         srcRedImg.setFilename(red_name);
00616         int imgRedNr = pano.addImage(srcRedImg);
00617         lenses.updatePartNumbers();
00618         lenses.switchParts(imgRedNr, 0);
00619 
00620         std::string green_name;
00621         if( g_param.green_name.size())
00622         {
00623             green_name=g_param.green_name;
00624         }
00625         else
00626         {
00627             green_name=std::string("green_")+filename;
00628         }
00629 
00630         HuginBase::SrcPanoImage srcGreenImg;
00631         srcGreenImg.setSize(imgInfo.size());
00632         srcGreenImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00633         srcGreenImg.setHFOV(10);
00634         srcGreenImg.setCropFactor(1);
00635         srcGreenImg.setFilename(green_name);
00636         int imgGreenNr = pano.addImage(srcGreenImg);
00637         lenses.updatePartNumbers();
00638         lenses.switchParts(imgGreenNr, 0);
00639 
00640         std::string blue_name;
00641         if( g_param.blue_name.size())
00642         {
00643             blue_name=g_param.blue_name;
00644         }
00645         else
00646         {
00647             blue_name=std::string("blue_")+filename;
00648         }
00649 
00650         HuginBase::SrcPanoImage srcBlueImg;
00651         srcBlueImg.setSize(imgInfo.size());
00652         srcBlueImg.setProjection(HuginBase::SrcPanoImage::RECTILINEAR);
00653         srcBlueImg.setHFOV(10);
00654         srcBlueImg.setCropFactor(1);
00655         srcBlueImg.setFilename(blue_name);
00656         int imgBlueNr = pano.addImage(srcBlueImg);
00657         lenses.updatePartNumbers();
00658         lenses.switchParts(imgBlueNr, 0);
00659 
00660         // lens variables are linked by default. Unlink the field of view and
00661         // the radial distortion.
00662         lenses.unlinkVariablePart(HuginBase::ImageVariableGroup::IVE_HFOV, 0);
00663         lenses.unlinkVariablePart(HuginBase::ImageVariableGroup::IVE_RadialDistortion, 0);
00664 
00665         // setup output to be exactly similar to input image
00666         HuginBase::PanoramaOptions opts;
00667         opts.setProjection(HuginBase::PanoramaOptions::RECTILINEAR);
00668         opts.setHFOV(srcGreenImg.getHFOV(), false);
00669         opts.setWidth(srcGreenImg.getSize().x, false);
00670         opts.setHeight(srcGreenImg.getSize().y);
00671 
00672         // output to tiff format
00673         opts.outputFormat = HuginBase::PanoramaOptions::TIFF_m;
00674         opts.tiff_saveROI = false;
00675         // m estimator, to be more robust against points on moving objects
00676         opts.huberSigma = 0.5;
00677         pano.setOptions(opts);
00678 
00679         createCtrlPoints(pano, imgOrig, imgRedNr, imgGreenNr, imgBlueNr, g_param.scale, g_param.nPoints, g_param.grid);
00680 
00681         main2(pano);
00682     }
00683     catch (std::exception& e)
00684     {
00685         std::cerr << "ERROR: caught exception: " << e.what() << std::endl;
00686         return 1;
00687     }
00688     return 0;
00689 }
00690 
00691 int processPTO(const char* filename)
00692 {
00693     HuginBase::Panorama pano;
00694 
00695     std::ifstream ptofile(filename);
00696     if (ptofile.bad())
00697     {
00698         std::cerr << "could not open script : " << filename << std::endl;
00699         return 1;
00700     }
00701     pano.setFilePrefix(hugin_utils::getPathPrefix(filename));
00702     AppBase::DocumentData::ReadWriteError err = pano.readData(ptofile);
00703     if (err != AppBase::DocumentData::SUCCESSFUL)
00704     {
00705         std::cerr << "error while parsing script: " << filename << std::endl;
00706         return 1;
00707     }
00708 
00709     return main2(pano);
00710 }
00711 
00712 void resetValues(HuginBase::Panorama& pano)
00713 {
00714     for (unsigned int i=0; i < 3; i++)
00715     {
00716         HuginBase::SrcPanoImage img = pano.getSrcImage(i);
00717 
00718         img.setHFOV(10);
00719 
00720         std::vector<double> radialDist(4);
00721         radialDist[0] = 0;
00722         radialDist[1] = 0;
00723         radialDist[2] = 0;
00724         radialDist[3] = 1;
00725         img.setRadialDistortion(radialDist);
00726 
00727         img.setRadialDistortionCenterShift(hugin_utils::FDiff2D(0, 0));
00728 
00729         pano.setSrcImage(i, img);
00730     }
00731 }
00732 
00733 void print_result(HuginBase::Panorama& pano)
00734 {
00735     double dist[3][3]; // a,b,c for all imgs
00736     double shift[2];   // x,y shift
00737     double hfov[3];
00738 
00739     for (unsigned int i=0; i < 3; i++)
00740     {
00741         HuginBase::SrcPanoImage img = pano.getSrcImage(i);
00742         hfov[i] = img.getHFOV();
00743         dist[i][0] = img.getRadialDistortion()[0];
00744         dist[i][1] = img.getRadialDistortion()[1];
00745         dist[i][2] = img.getRadialDistortion()[2];
00746         if (i == 0)
00747         {
00748             shift[0] = img.getRadialDistortionCenterShift().x;
00749             shift[1] = img.getRadialDistortionCenterShift().y;
00750         }
00751     }
00752 
00753     /* compute new a,b,c,d from a,b,c,v */
00754     double distnew[3][4];
00755     for (unsigned int i=0 ; i<3 ; i++)
00756     {
00757         double scale = hfov[1] / hfov[i];
00758         for (unsigned int j=0 ; j<3 ; j++)
00759         {
00760             distnew[i][j] = dist[i][j]*pow(scale, (int)(4-j));
00761         }
00762         distnew[i][3] = scale*(1 - dist[i][0] - dist[i][1] - dist[i][2]);
00763     }
00764 
00765     if (hugin_utils::roundi(shift[0]) == 0 && hugin_utils::roundi(shift[1]) == 0)
00766     {
00767         fprintf(stdout, "-r %.7f:%.7f:%.7f:%.7f "
00768             "-b %.7f:%.7f:%.7f:%.7f ",
00769             distnew[0][0], distnew[0][1], distnew[0][2], distnew[0][3],
00770             distnew[2][0], distnew[2][1], distnew[2][2], distnew[2][3]);
00771     }
00772     else
00773     {
00774         fprintf(stdout, "-r %.7f:%.7f:%.7f:%.7f "
00775             "-b %.7f:%.7f:%.7f:%.7f "
00776             "-x %d:%d\n",
00777             distnew[0][0], distnew[0][1], distnew[0][2], distnew[0][3],
00778             distnew[2][0], distnew[2][1], distnew[2][2], distnew[2][3],
00779             hugin_utils::roundi(shift[0]), hugin_utils::roundi(shift[1]));
00780     };
00781     if (g_param.saveDB)
00782     {
00783         // save information into database
00784         // read EXIF data, using routines of HuginBase::SrcPanoImage
00785         HuginBase::SrcPanoImage img;
00786         img.setFilename(g_param.basename);
00787         img.readEXIF();
00788         std::string lensname = img.getDBLensName();
00789         if (lensname.empty())
00790         {
00791             // no suitable lensname found in exif data, ask user
00792             std::cout << std::endl << "For saving information tca data into database no suitable information" << std::endl
00793                 << "found in EXIF data." << std::endl
00794                 << "For fixed lens cameras leave lensname empty." << std::endl << std::endl << "Lensname: ";
00795             char input[256];
00796             std::cin.getline(input, 255);
00797             lensname = hugin_utils::StrTrim(std::string(input));
00798             if (lensname.empty())
00799             {
00800                 std::string camMaker;
00801                 std::string camModel;
00802                 std::cout << "Camera maker: ";
00803                 while (camMaker.empty())
00804                 {
00805                     std::cin.getline(input, 255);
00806                     camMaker = hugin_utils::StrTrim(std::string(input));
00807                 };
00808                 std::cout << "Camera model: ";
00809                 while (camModel.empty())
00810                 {
00811                     std::cin.getline(input, 255);
00812                     camModel = hugin_utils::StrTrim(std::string(input));
00813                 };
00814                 lensname = camMaker.append("|").append(camModel);
00815             };
00816         };
00817         double focal = img.getExifFocalLength();
00818         if (fabs(focal) < 0.1f)
00819         {
00820             std::cout << "Real focal length (in mm): ";
00821             while (fabs(focal) < 0.1f)
00822             {
00823                 std::cin >> focal;
00824                 focal = fabs(focal);
00825             };
00826         };
00827         HuginBase::LensDB::LensDB& lensDB = HuginBase::LensDB::LensDB::GetSingleton();
00828         std::vector<double> redDistData(4, 0.0f);
00829         std::vector<double> blueDistData(4, 0.0f);
00830         redDistData[0] = distnew[0][0];
00831         redDistData[1] = distnew[0][1];
00832         redDistData[2] = distnew[0][2];
00833         redDistData[3] = distnew[0][3];
00834         blueDistData[0] = distnew[2][0];
00835         blueDistData[1] = distnew[2][1];
00836         blueDistData[2] = distnew[2][2];
00837         blueDistData[3] = distnew[2][3];
00838         if (lensDB.SaveTCA(lensname, focal, redDistData, blueDistData))
00839         {
00840             std::cout << std::endl << std::endl << "TCA data for " << lensname << " @ " << focal << " mm successful saved into database." << std::endl;
00841         }
00842         else
00843         {
00844             std::cout << std::endl << std::endl << "Could not save data into database." << std::endl;
00845         };
00846         HuginBase::LensDB::LensDB::Clean();
00847     }
00848 }
00849 
00850 int main2(HuginBase::Panorama& pano)
00851 {
00852     if (g_param.reset)
00853     {
00854         resetValues(pano);
00855     }
00856 
00857     for (int i = 0; i < 10; i++)
00858     {
00859         if (g_param.optMethod == 0)
00860         {
00861             optimize_old(pano);
00862         }
00863         else if (g_param.optMethod == 1)
00864         {
00865             optimize_new(pano);
00866         }
00867         else
00868         {
00869             std::cerr << "Unknown optimizer strategy." << std::endl
00870                 << "Using newfit method." << std::endl;
00871             g_param.optMethod = 1;
00872             optimize_new(pano);
00873         };
00874 
00875         HuginBase::CPVector cps = pano.getCtrlPoints();
00876         HuginBase::CPVector newCPs;
00877         for (HuginBase::CPVector::const_iterator it = cps.begin(); it != cps.end(); ++it)
00878         {
00879             if (it->error < g_param.cpErrorThreshold)
00880             {
00881                 newCPs.push_back(*it);
00882             }
00883         }
00884         if (g_param.verbose > 0)
00885         {
00886             std::cerr << "Ctrl points before pruning: " << cps.size() << ", after: " << newCPs.size() << std::endl;
00887         }
00888         pano.setCtrlPoints(newCPs);
00889 
00890         if (cps.size() == newCPs.size())
00891         {
00892             // no points were removed, do not re-optimize
00893             break;
00894         }
00895     }
00896 
00897     if (!g_param.ptoOutputFile.empty())
00898     {
00899         HuginBase::OptimizeVector optvars;
00900         get_optvars(optvars);
00901         HuginBase::UIntSet allImgs;
00902         fill_set(allImgs, 0, pano.getNrOfImages() - 1);
00903         std::ofstream script(g_param.ptoOutputFile.c_str());
00904         pano.printPanoramaScript(script, optvars, pano.getOptions(), allImgs, true, "");
00905     }
00906 
00907     print_result(pano);
00908     return 0;
00909 }
00910 
00911 int main(int argc, char* argv[])
00912 {
00913     // parse arguments
00914     const char* optstring = "hlm:o:rt:vw:R:G:B:s:g:n:";
00915     enum
00916     {
00917         SAVEDATABASE=1000
00918     };
00919     static struct option longOptions[] =
00920     {
00921         { "save-into-database", no_argument, NULL, SAVEDATABASE },
00922         { "help", no_argument, NULL, 'h' },
00923         0
00924     };
00925     int c;
00926     bool parameter_request_seen=false;
00927     int optionIndex = 0;
00928     opterr = 0;
00929 
00930     while ((c = getopt_long(argc, argv, optstring, longOptions, &optionIndex)) != -1)
00931     {
00932         switch (c)
00933         {
00934             case 'h':
00935                 usage(hugin_utils::stripPath(argv[0]).c_str());
00936                 return 0;
00937             case 'l':
00938                 g_param.load = true;
00939                 break;
00940             case 'm':
00941                 g_param.optMethod = atoi(optarg);
00942                 break;
00943             case 'o':
00944                 {
00945                     char* optptr = optarg;
00946                     while (*optptr != 0)
00947                     {
00948                         if ((*optptr == 'a') || (*optptr == 'b') ||
00949                                 (*optptr == 'c') || (*optptr == 'v') ||
00950                                 (*optptr == 'd') || (*optptr == 'e'))
00951                         {
00952                             g_param.optvars.insert(std::string(optptr, 1));
00953                         }
00954                         optptr++;
00955                     }
00956                     parameter_request_seen=true;
00957                 }
00958                 break;
00959             case 'r':
00960                 g_param.reset = true;
00961                 break;
00962             case 't':
00963                 g_param.cpErrorThreshold = atof(optarg);
00964                 if (g_param.cpErrorThreshold <= 0)
00965                 {
00966                     std::cerr << "Invalid parameter: control point error threshold (-t) must be greater than 0" << std::endl;
00967                     return 1;
00968                 }
00969                 break;
00970             case 'v':
00971                 g_param.verbose++;
00972                 break;
00973             case 'w':
00974                 g_param.ptoOutputFile = optarg;
00975                 break;
00976             case 'R':
00977                 g_param.red_name = optarg;
00978                 break;
00979             case 'G':
00980                 g_param.green_name = optarg;
00981                 break;
00982             case 'B':
00983                 g_param.blue_name = optarg;
00984                 break;
00985             case 's':
00986                 g_param.scale=atof( optarg);
00987                 break;
00988             case 'n':
00989                 g_param.nPoints=atoi( optarg);
00990                 break;
00991             case 'g':
00992                 g_param.grid=atoi(optarg);
00993                 break;
00994             case SAVEDATABASE:
00995                 g_param.saveDB = true;
00996                 break;
00997             default:
00998                 std::cerr << "Invalid parameter: '" << argv[optind-1] << " " << optarg << "'" << std::endl;
00999                 usage(hugin_utils::stripPath(argv[0]).c_str());
01000                 return 1;
01001         }
01002     };
01003 
01004     if ((argc - optind) != 1)
01005     {
01006         usage(hugin_utils::stripPath(argv[0]).c_str());
01007         return 1;
01008     }
01009 
01010     // If no parameters were requested to be optimised, we optimize the
01011     // default parameters.
01012     if ( !parameter_request_seen)
01013     {
01014         for ( const char* dop=DEFAULT_OPTIMISATION_PARAMETER;
01015                 *dop != 0; ++dop)
01016         {
01017             g_param.optvars.insert( std::string( dop, 1));
01018         }
01019     }
01020 
01021     // Program will crash if nothing is to be optimised.
01022     if ( g_param.optvars.empty())
01023     {
01024         std::cerr << "No parameters to optimize." << std::endl;
01025         usage(hugin_utils::stripPath(argv[0]).c_str());
01026         return 1;
01027     }
01028 
01029     if (!g_param.load)
01030     {
01031         g_param.basename = argv[optind];
01032         vigra::ImageImportInfo firstImgInfo(argv[optind]);
01033         std::string pixelType = firstImgInfo.getPixelType();
01034         if (pixelType == "UINT8")
01035         {
01036             return processImg<vigra::RGBValue<vigra::UInt8> >(argv[optind]);
01037         }
01038         else if (pixelType == "INT16")
01039         {
01040             return processImg<vigra::RGBValue<vigra::Int16> >(argv[optind]);
01041         }
01042         else if (pixelType == "UINT16")
01043         {
01044             return processImg<vigra::RGBValue<vigra::UInt16> >(argv[optind]);
01045         }
01046         else
01047         {
01048             std::cerr << " ERROR: unsupported pixel type: " << pixelType << std::endl;
01049             return 1;
01050         }
01051     }
01052     else
01053     {
01054         return processPTO(argv[optind]);
01055     }
01056 
01057     return 0;
01058 }
01059 

Generated on 28 May 2016 for Hugintrunk by  doxygen 1.4.7