FindLines.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00008 /***************************************************************************
00009  *   Copyright (C) 2009 by Tim Nugent                                      *
00010  *   timnugent@gmail.com                                                   *
00011  *                                                                         *
00012  *   This program is free software; you can redistribute it and/or modify  *
00013  *   it under the terms of the GNU General Public License as published by  *
00014  *   the Free Software Foundation; either version 2 of the License, or     *
00015  *   (at your option) any later version.                                   *
00016  *                                                                         *
00017  *   This program is distributed in the hope that it will be useful,       *
00018  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00019  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00020  *   GNU General Public License for more details.                          *
00021  *                                                                         *
00022  *   You should have received a copy of the GNU General Public License     *
00023  *   along with this program.  If not, see <http://www.gnu.org/licenses/>. *
00024  ***************************************************************************/
00025 
00026 #include "vigra/edgedetection.hxx"
00027 #include "FindLines.h"
00028 #include "FindN8Lines.h"
00029 #include <algorithms/nona/FitPanorama.h>
00030 #include <algorithms/basic/CalculateOptimalROI.h>
00031 #include <nona/RemappedPanoImage.h>
00032 #include <algorithms/optimizer/PTOptimizer.h>
00033 #include "algorithms/basic/CalculateCPStatistics.h"
00034 
00035 namespace HuginLines
00036 {
00037 
00038 template <class SrcImageIterator, class SrcAccessor, class DestImage>
00039 double resize_image(const vigra::triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src, DestImage& dest, int resize_dimension)
00040 {
00041     // Re-size to max dimension
00042     double sizefactor=1.0;
00043     const vigra::Size2D inputSize(src.second - src.first);
00044     if (inputSize.width() > resize_dimension || inputSize.height() > resize_dimension)
00045     {
00046         int nw;
00047         int nh;
00048         if (inputSize.width() >= inputSize.height())
00049         {
00050             sizefactor = (double)resize_dimension / inputSize.width();
00051             // calculate new image size
00052             nw = resize_dimension;
00053             nh = static_cast<int>(0.5 + (sizefactor*inputSize.height()));
00054         }
00055         else
00056         {
00057             sizefactor = (double)resize_dimension / inputSize.height();
00058             // calculate new image size
00059             nw = static_cast<int>(0.5 + (sizefactor*inputSize.width()));
00060             nh = resize_dimension;
00061         }
00062 
00063         // create an image of appropriate size
00064         dest.resize(nw, nh);
00065         // resize the image, using a bi-cubic spline algorithm
00066         resizeImageNoInterpolation(src, destImageRange(dest));
00067     }
00068     else
00069     {
00070         dest.resize(inputSize);
00071         copyImage(src, destImage(dest));
00072     };
00073     return 1.0/sizefactor;
00074 }
00075 
00076 vigra::BImage* detectEdges(const vigra::UInt8RGBImage& input, const double scale, const double threshold, const unsigned int resize_dimension, double& size_factor)
00077 {
00078     // Resize image
00079     vigra::UInt8Image scaled;
00080     size_factor = resize_image(vigra::srcImageRange(input, vigra::RGBToGrayAccessor<vigra::RGBValue<vigra::UInt8> >()), scaled, resize_dimension);
00081 
00082     // Run Canny edge detector
00083     vigra::BImage* image = new vigra::BImage(scaled.width(), scaled.height(), 255);
00084     vigra::cannyEdgeImage(vigra::srcImageRange(scaled), vigra::destImage(*image), scale, threshold, 0);
00085     return image;
00086 };
00087 
00088 vigra::BImage* detectEdges(const vigra::BImage& input, const double scale, const double threshold, const unsigned int resize_dimension, double& size_factor)
00089 {
00090     // Resize image
00091     vigra::UInt8Image scaled;
00092     size_factor=resize_image(vigra::srcImageRange(input), scaled, resize_dimension);
00093 
00094     // Run Canny edge detector
00095     vigra::BImage* image = new vigra::BImage(scaled.width(), scaled.height(), 255);
00096     vigra::cannyEdgeImage(vigra::srcImageRange(scaled), vigra::destImage(*image), scale, threshold, 0);
00097     return image;
00098 };
00099 
00100 double calculate_focal_length_pixels(double focal_length,double cropFactor,double width, double height)
00101 {
00102     double pixels_per_mm = 0;
00103     if (cropFactor > 1)
00104     {
00105         pixels_per_mm= (cropFactor/24.0)* ((width>height)?height:width);
00106     }
00107     else
00108     {
00109         pixels_per_mm= (24.0/cropFactor)* ((width>height)?height:width);
00110     }
00111     return focal_length*pixels_per_mm;
00112 }
00113 
00114 
00115 Lines findLines(vigra::BImage& edge, double length_threshold, double focal_length,double crop_factor)
00116 {
00117     unsigned int longest_dimension=(edge.width() > edge.height()) ? edge.width() : edge.height();
00118     double min_line_length_squared=(length_threshold*longest_dimension)*(length_threshold*longest_dimension);
00119 
00120     int lmin = int(sqrt(min_line_length_squared));
00121     double flpix=calculate_focal_length_pixels(focal_length,crop_factor,edge.width(),edge.height());
00122 
00123     vigra::BImage lineImage = edgeMap2linePts(edge);
00124     Lines lines;
00125     linePts2lineList( lineImage, lmin, flpix, lines );
00126 
00127     return lines;
00128 };
00129 
00130 void ScaleLines(Lines& lines,const double scale)
00131 {
00132     for(unsigned int i=0; i<lines.size(); i++)
00133     {
00134         for(unsigned int j=0; j<lines[i].line.size(); j++)
00135         {
00136             lines[i].line[j]*=scale;
00137         };
00138     };
00139 };
00140 
00141 HuginBase::CPVector GetControlPoints(const SingleLine& line,const unsigned int imgNr, const unsigned int lineNr,const unsigned int numberOfCtrlPoints)
00142 {
00143     HuginBase::CPVector cpv;
00144     double interval = (line.line.size()-1)/(1.0*numberOfCtrlPoints);
00145     for(unsigned int k = 0; k < numberOfCtrlPoints; k++)
00146     {
00147         int start = (int)(k * interval);
00148         int stop =  (int)((k+1) * interval);
00149         HuginBase::ControlPoint cp(imgNr,line.line[start].x, line.line[start].y,
00150                                    imgNr,line.line[stop].x, line.line[stop].y,lineNr);
00151         cpv.push_back(cp);
00152     };
00153     return cpv;
00154 };
00155 
00156 #define MAX_RESIZE_DIM 1600
00157 
00158 struct VerticalLine
00159 {
00160     vigra::Point2D start;
00161     vigra::Point2D end;
00162 };
00163 
00164 typedef std::vector<VerticalLine> VerticalLineVector;
00165 
00166 //return footpoint of point p on line between point p1 and p2
00167 vigra::Point2D GetFootpoint(vigra::Point2D p, vigra::Point2D p1, vigra::Point2D p2)
00168 {
00169     hugin_utils::FDiff2D diff=p2-p1;
00170     double u=((p.x-p1.x)*(p2.x-p1.x)+(p.y-p1.y)*(p2.y-p1.y))/hugin_utils::sqr(hugin_utils::norm(diff));
00171     diff*=u;
00172     return vigra::Point2D(p1.x+diff.x,p1.y+diff.y);
00173 };
00174 
00175 //linear fit of given line, returns endpoints of fitted line
00176 VerticalLine FitLine(SingleLine line)
00177 {
00178     size_t n=line.line.size();
00179     VerticalLine vl;
00180     double s_x=0;
00181     double s_y=0;
00182     double s_xy=0;
00183     double s_x2=0;
00184     for(size_t i=0;i<n;i++)
00185     {
00186         s_x+=(double)line.line[i].x/n;
00187         s_y+=(double)line.line[i].y/n;
00188         s_xy+=(double)line.line[i].x*line.line[i].y/n;
00189         s_x2+=(double)line.line[i].x*line.line[i].x/n;
00190     };
00191     if(std::abs(s_x2-s_x*s_x)<0.00001)
00192     {
00193         //vertical line needs special treatment
00194         vl.start.x=s_x;
00195         vl.start.y=line.line[0].y;
00196         vl.end.x=s_x;
00197         vl.end.y=line.line[n-1].y;
00198     }
00199     else
00200     {
00201         //calculate slope and offset
00202         double slope=(s_xy-s_x*s_y)/(s_x2-s_x*s_x);
00203         double offset=s_y-slope*s_x;
00204         //convert to parametric form
00205         vigra::Point2D p1(0,offset);
00206         vigra::Point2D p2(100,100*slope+offset);
00207         //calculate footpoint of first and last point
00208         vl.start=GetFootpoint(line.line[0],p1,p2);
00209         vl.end=GetFootpoint(line.line[n-1],p1,p2);
00210     };
00211     return vl;
00212 };
00213 
00214 //filter detected lines
00215 //return fitted lines which have only a small deviation from the vertical
00216 VerticalLineVector FilterLines(Lines lines,double roll)
00217 {
00218     VerticalLineVector vertLines;
00219     if(lines.size()>0)
00220     {
00221         for(Lines::const_iterator it=lines.begin(); it!=lines.end(); ++it)
00222         {
00223             if((*it).status==valid_line && (*it).line.size()>2)
00224             {
00225                 VerticalLine vl=FitLine(*it);
00226                 vigra::Diff2D diff=vl.end-vl.start;
00227                 if(diff.magnitude()>20)
00228                 {
00229                     if(std::abs((diff.x*cos(DEG_TO_RAD(roll))+diff.y*sin(DEG_TO_RAD(roll)))/diff.magnitude())<0.1)
00230                     {
00231                         vertLines.push_back(vl);
00232                     };
00233                 };
00234             };
00235         };
00236     };
00237     return vertLines;
00238 };
00239 
00240 //function to sort HuginBase::CPVector by error distance
00241 bool SortByError(const HuginBase::ControlPoint& cp1, const HuginBase::ControlPoint& cp2)
00242 {
00243     return cp1.error<cp2.error;
00244 };
00245 
00246 template <class ImageType>
00247 HuginBase::CPVector _getVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,ImageType& image, const unsigned int nrLines)
00248 {
00249     HuginBase::CPVector verticalLines;
00250     HuginBase::CPVector detectedLines;
00251     const HuginBase::SrcPanoImage& srcImage=pano.getImage(imgNr);
00252     bool needsRemap=srcImage.getProjection()!=HuginBase::SrcPanoImage::RECTILINEAR;
00253     double roll=(needsRemap?0:srcImage.getRoll());
00254     double size_factor=1.0;
00255     HuginBase::SrcPanoImage remappedImage;
00256     HuginBase::PanoramaOptions opts;
00257     vigra::BImage* edge;
00258     if(!needsRemap)
00259     {
00260         //rectilinear image can be used as is
00261         //detect edges
00262         edge=detectEdges(image,2,4,MAX_RESIZE_DIM,size_factor);
00263     }
00264     else
00265     {
00266         //remap all other image to equirectangular
00267         //create temporary SrcPanoImage, set appropriate image variables
00268         remappedImage=pano.getSrcImage(imgNr);
00269         remappedImage.setYaw(0);
00270         remappedImage.setPitch(0);
00271         remappedImage.setX(0);
00272         remappedImage.setY(0);
00273         remappedImage.setZ(0);
00274         remappedImage.setExposureValue(0);
00275         remappedImage.setEMoRParams(std::vector<float>(5, 0.0));
00276         remappedImage.deleteAllMasks();
00277         remappedImage.setActive(true);
00278         //create PanoramaOptions for remapping of image
00279         opts.setProjection(HuginBase::PanoramaOptions::EQUIRECTANGULAR);
00280         opts.setWidth(MAX_RESIZE_DIM);
00281         opts.outputExposureValue=0;
00282         //calculate output canvas size
00283         HuginBase::Panorama tempPano;
00284         tempPano.addImage(remappedImage);
00285         tempPano.setOptions(opts);
00286 
00287         HuginBase::CalculateFitPanorama fitPano(tempPano);
00288         fitPano.run();
00289         opts.setHFOV(fitPano.getResultHorizontalFOV());
00290         opts.setHeight(hugin_utils::roundi(fitPano.getResultHeight()));
00291         tempPano.setOptions(opts);
00292 
00293         //finally remap image
00294         HuginBase::Nona::RemappedPanoImage<ImageType,vigra::BImage>* remapped=new HuginBase::Nona::RemappedPanoImage<ImageType,vigra::BImage>;
00295         AppBase::ProgressDisplay* progress=new AppBase::DummyProgressDisplay();
00296         remapped->setPanoImage(remappedImage,opts,opts.getROI());
00297         remapped->remapImage(vigra::srcImageRange(image),vigra_ext::INTERP_CUBIC, progress);
00298         ImageType remappedBitmap=remapped->m_image;
00299         //detect edges
00300         edge=detectEdges(remappedBitmap,2,4,std::max(remappedBitmap.width(),remappedBitmap.height())+10,size_factor);
00301         delete remapped;
00302         delete progress;
00303     };
00304     //detect lines
00305     //we need the focal length
00306     double focalLength=srcImage.getExifFocalLength();
00307     if(focalLength==0)
00308     {
00309         focalLength=HuginBase::SrcPanoImage::calcFocalLength(
00310             srcImage.getProjection(),srcImage.getHFOV(),srcImage.getCropFactor(),srcImage.getSize());
00311     };
00312     Lines foundLines=findLines(*edge,0.05,focalLength,srcImage.getCropFactor());
00313     delete edge;
00314     //filter results
00315     VerticalLineVector filteredLines=FilterLines(foundLines,roll);
00316     //create control points
00317     if(filteredLines.size()>0)
00318     {
00319         //we need to transform the coordinates to image coordinates because the detection
00320         //worked on smaller images or in remapped image
00321         HuginBase::PTools::Transform transform;
00322         if(needsRemap)
00323         {
00324             transform.createTransform(remappedImage,opts);
00325         };
00326         for(size_t i=0; i<filteredLines.size(); i++)
00327         {
00328             HuginBase::ControlPoint cp;
00329             cp.image1Nr=0;
00330             cp.image2Nr=0;
00331             cp.mode=HuginBase::ControlPoint::X;
00332             if(!needsRemap)
00333             {
00334                 cp.x1=filteredLines[i].start.x*size_factor;
00335                 cp.y1=filteredLines[i].start.y*size_factor;
00336                 cp.x2=filteredLines[i].end.x*size_factor;
00337                 cp.y2=filteredLines[i].end.y*size_factor;
00338             }
00339             else
00340             {
00341                 double xout;
00342                 double yout;
00343                 if(!transform.transformImgCoord(xout,yout,filteredLines[i].start.x,filteredLines[i].start.y))
00344                 {
00345                     continue;
00346                 };
00347                 cp.x1=xout;
00348                 cp.y1=yout;
00349                 if(!transform.transformImgCoord(xout,yout,filteredLines[i].end.x,filteredLines[i].end.y))
00350                 {
00351                     continue;
00352                 };
00353                 cp.x2=xout;
00354                 cp.y2=yout;
00355             };
00356             if(cp.x1>=0 && cp.x1<srcImage.getWidth() && cp.y1>=0 && cp.y1<srcImage.getHeight() &&
00357                cp.x2>=0 && cp.x2<srcImage.getWidth() && cp.y2>=0 && cp.y2<srcImage.getHeight())
00358             {
00359                 detectedLines.push_back(cp);
00360             };
00361         };
00362         //now a final check of the found vertical lines
00363         //we optimize the pano with a single image and disregard vertical lines with bigger errors
00364         //we need at least 2 lines
00365         if(detectedLines.size()>1)
00366         {
00367             HuginBase::Panorama tempPano;
00368             HuginBase::SrcPanoImage tempImage=pano.getSrcImage(imgNr);
00369             tempImage.setYaw(0);
00370             tempImage.setPitch(0);
00371             tempImage.setRoll(0);
00372             tempImage.setX(0);
00373             tempImage.setY(0);
00374             tempImage.setZ(0);
00375             tempPano.addImage(tempImage);
00376             for(size_t i=0; i<detectedLines.size(); i++)
00377             {
00378                 tempPano.addCtrlPoint(detectedLines[i]);
00379             };
00380             HuginBase::PanoramaOptions opt2;
00381             opt2.setProjection(HuginBase::PanoramaOptions::EQUIRECTANGULAR);
00382             tempPano.setOptions(opt2);
00383             HuginBase::OptimizeVector optVec;
00384             std::set<std::string> imgopt;
00385             imgopt.insert("p");
00386             imgopt.insert("r");
00387             optVec.push_back(imgopt);
00388             tempPano.setOptimizeVector(optVec);
00389             // ARGH the panotools optimizer uses global variables is not reentrant
00390 #pragma omp critical
00391             {
00392                 HuginBase::PTools::optimize(tempPano);
00393             }
00394             //first filter stage
00395             //we disregard all lines with big error
00396             //calculate statistic and determine limit
00397             double minError,maxError,mean,var;
00398             HuginBase::CalculateCPStatisticsError::calcCtrlPntsErrorStats(tempPano,minError,maxError,mean,var);
00399             detectedLines=tempPano.getCtrlPoints();
00400             double limit=mean+sqrt(var);
00401             maxError=0;
00402             for(int i=detectedLines.size()-1; i>=0; i--)
00403             {
00404                 if(detectedLines[i].error>limit)
00405                 {
00406                     detectedLines.erase(detectedLines.begin()+i);
00407                 }
00408                 else
00409                 {
00410                     //we need the max error of the remaining lines for the next step
00411                     maxError=std::max(detectedLines[i].error,maxError);
00412                 };
00413             };
00414             if(detectedLines.size()>0 && maxError>0) //security check, should never be false
00415             {
00416                 //now keep only the best nrLines lines
00417                 //we are using error and line length as figure of merrit
00418                 for(size_t i=0;i<detectedLines.size();i++)
00419                 {
00420                     double length=sqrt(hugin_utils::sqr(detectedLines[i].x2-detectedLines[i].x1)+hugin_utils::sqr(detectedLines[i].y2-detectedLines[i].y1));
00421                     //calculate number of merrit
00422                     detectedLines[i].error=detectedLines[i].error/maxError+(1.0-std::min(length,500.0)/500.0);
00423                 };
00424                 std::sort(detectedLines.begin(),detectedLines.end(),SortByError);
00425                 //only save best nrLines control points
00426                 for(size_t i=0;i<detectedLines.size() && i<nrLines; i++)
00427                 {
00428                     HuginBase::ControlPoint cp=detectedLines[i];
00429                     cp.image1Nr=imgNr;
00430                     cp.image2Nr=imgNr;
00431                     cp.error=0;
00432                     verticalLines.push_back(cp);
00433                 };
00434             };
00435         }
00436         else
00437         {
00438             //if only one line was detected we do a special check
00439             //the allow deviation between line and roll angle is checked more narrow than in the first check
00440             if(detectedLines.size()==1)
00441             {
00442                 vigra::Diff2D diff((double)detectedLines[0].x2-detectedLines[0].x1,(double)detectedLines[0].y2-detectedLines[0].y1);
00443                 if(std::abs((diff.x*cos(DEG_TO_RAD(roll))+diff.y*sin(DEG_TO_RAD(roll)))/diff.magnitude())<0.05)
00444                 {
00445                     HuginBase::ControlPoint cp=detectedLines[0];
00446                     cp.image1Nr=imgNr;
00447                     cp.image2Nr=imgNr;
00448                     cp.error=0;
00449                     verticalLines.push_back(cp);
00450                 };
00451             };
00452         };
00453     };
00454     return verticalLines;
00455 };
00456 
00457 HuginBase::CPVector GetVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,vigra::UInt8RGBImage& image, const unsigned int nrLines)
00458 {
00459     return _getVerticalLines(pano, imgNr, image, nrLines);
00460 };
00461 
00462 HuginBase::CPVector GetVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,vigra::BImage& image, const unsigned int nrLines)
00463 {
00464     return _getVerticalLines(pano, imgNr, image, nrLines);
00465 };
00466 
00467 }; //namespace

Generated on 1 Jul 2016 for Hugintrunk by  doxygen 1.4.7