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

Generated on Tue Sep 16 01:25:40 2014 for Hugintrunk by  doxygen 1.3.9.1