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

Generated on Tue Sep 2 01:25:42 2014 for Hugintrunk by  doxygen 1.3.9.1