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

Generated on 31 Aug 2015 for Hugintrunk by  doxygen 1.4.7