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 //return footpoint of point p on line between point p1 and p2
00159 vigra::Point2D GetFootpoint(const vigra::Point2D& p, const vigra::Point2D& p1, const vigra::Point2D& p2, double& u)
00160 {
00161     hugin_utils::FDiff2D diff = p2 - p1;
00162     u = ((p.x - p1.x)*(p2.x - p1.x) + (p.y - p1.y)*(p2.y - p1.y)) / hugin_utils::sqr(hugin_utils::norm(diff));
00163     diff *= u;
00164     return vigra::Point2D(p1.x + diff.x, p1.y + diff.y);
00165 };
00166 
00167 vigra::Point2D GetFootpoint(const vigra::Point2D& p, const vigra::Point2D& p1, const vigra::Point2D& p2)
00168 {
00169     double u;
00170     return GetFootpoint(p, p1, p2, u);
00171 };
00172 
00173 class VerticalLine
00174 {
00175 public:
00176     void SetStart(const vigra::Point2D point)
00177     {
00178         m_start = point;
00179     };
00180     void SetStart(const int x, const int y)
00181     {
00182         m_start = vigra::Point2D(x, y);
00183     };
00184     void SetEnd(const vigra::Point2D point)
00185     {
00186         m_end = point;
00187     };
00188     void SetEnd(const int x, const int y)
00189     {
00190         m_end = vigra::Point2D(x, y);
00191     };
00192     double GetLineLength() const
00193     {
00194         return (m_end - m_start).magnitude();
00195     };
00196     double GetEstimatedDistance(const VerticalLine& otherLine) const
00197     {
00198         auto getDist = [](const vigra::Point2D& p, const vigra::Point2D& p1, const vigra::Point2D& p2)->double
00199         {
00200             double t;
00201             const vigra::Point2D endPoint = GetFootpoint(p, p1, p2, t);
00202             if (-0.1 < t && t < 1.1)
00203             {
00204                 return (endPoint - p).magnitude();
00205             }
00206             else
00207             {
00208                 return DBL_MAX;
00209             };
00210         };
00211         return std::min({getDist(otherLine.GetStart(), m_start, m_end), getDist(otherLine.GetEnd(), m_start, m_end),
00212             getDist(m_start, otherLine.GetStart(), otherLine.GetEnd()), getDist(m_end, otherLine.GetStart(), otherLine.GetEnd())});
00213     }
00214     double GetAngle() const
00215     {
00216         return atan2(m_end.y - m_start.y, m_end.x - m_start.x);
00217     };
00218     const vigra::Point2D& GetStart() const
00219     {
00220         return m_start;
00221     };
00222     const vigra::Point2D& GetEnd() const
00223     {
00224         return m_end;
00225     };
00226 private:
00227     vigra::Point2D m_start;
00228     vigra::Point2D m_end;
00229 };
00230 
00231 typedef std::vector<VerticalLine> VerticalLineVector;
00232 
00233 //linear fit of given line, returns endpoints of fitted line
00234 VerticalLine FitLine(SingleLine line)
00235 {
00236     size_t n=line.line.size();
00237     VerticalLine vl;
00238     double s_x=0;
00239     double s_y=0;
00240     double s_xy=0;
00241     double s_x2=0;
00242     for(size_t i=0;i<n;i++)
00243     {
00244         s_x+=(double)line.line[i].x/n;
00245         s_y+=(double)line.line[i].y/n;
00246         s_xy+=(double)line.line[i].x*line.line[i].y/n;
00247         s_x2+=(double)line.line[i].x*line.line[i].x/n;
00248     };
00249     if(std::abs(s_x2-s_x*s_x)<0.00001)
00250     {
00251         //vertical line needs special treatment
00252         vl.SetStart(s_x, line.line[0].y);
00253         vl.SetEnd(s_x, line.line[n - 1].y);
00254     }
00255     else
00256     {
00257         //calculate slope and offset
00258         double slope=(s_xy-s_x*s_y)/(s_x2-s_x*s_x);
00259         double offset=s_y-slope*s_x;
00260         //convert to parametric form
00261         vigra::Point2D p1(0,offset);
00262         vigra::Point2D p2(100,100*slope+offset);
00263         //calculate footpoint of first and last point
00264         vl.SetStart(GetFootpoint(line.line[0], p1, p2));
00265         vl.SetEnd(GetFootpoint(line.line[n - 1], p1, p2));
00266     };
00267     return vl;
00268 };
00269 
00270 //filter detected lines
00271 //return fitted lines which have only a small deviation from the vertical
00272 VerticalLineVector FilterLines(Lines lines,double roll)
00273 {
00274     VerticalLineVector vertLines;
00275     if(!lines.empty())
00276     {
00277         for(Lines::const_iterator it=lines.begin(); it!=lines.end(); ++it)
00278         {
00279             if((*it).status==valid_line && (*it).line.size()>2)
00280             {
00281                 VerticalLine vl=FitLine(*it);
00282                 const vigra::Diff2D diff = vl.GetEnd() - vl.GetStart();
00283                 // check that line is long enough
00284                 if(diff.magnitude()>20)
00285                 {
00286                     // now check angle with respect to roll angle, accept only deviation of 5 (sin 5=0.1)
00287                     if(std::abs((diff.x*cos(DEG_TO_RAD(roll))+diff.y*sin(DEG_TO_RAD(roll)))/diff.magnitude())<0.1)
00288                     {
00289                         // check distance and angle to other lines
00290                         bool distanceBig = true;
00291                         for (auto& otherLine : vertLines)
00292                         {
00293                             // distance should be at least 80 pixel = 5 % from image width
00294                             if (vl.GetEstimatedDistance(otherLine) < 80)
00295                             {
00296                                 // now check if line are parallel = have the same angle (tan(3)=0.05)
00297                                 if (std::abs(vl.GetAngle() - otherLine.GetAngle()) < 0.05)
00298                                 {
00299                                     distanceBig = false;
00300                                     // both lines are close to each other, keep only the longer one
00301                                     if (vl.GetLineLength() > otherLine.GetLineLength())
00302                                     {
00303                                         otherLine = vl;
00304                                     }
00305                                     continue;
00306                                 };
00307                             };
00308                         }
00309                         if (distanceBig)
00310                         {
00311                             vertLines.push_back(vl);
00312                         };
00313                     };
00314                 };
00315             };
00316         };
00317     };
00318     return vertLines;
00319 };
00320 
00321 //function to sort HuginBase::CPVector by error distance
00322 bool SortByError(const HuginBase::ControlPoint& cp1, const HuginBase::ControlPoint& cp2)
00323 {
00324     return cp1.error<cp2.error;
00325 };
00326 
00327 class InvertedMaskAccessor
00328 {
00329 public:
00330     typedef vigra::UInt8 value_type;
00331     template <class ITERATOR>
00332     value_type operator()(ITERATOR const & i) const
00333     {
00334         return 255 - (*i);
00335     }
00336     template <class ITERATOR, class DIFFERENCE>
00337     value_type operator()(ITERATOR const & i, DIFFERENCE d) const
00338     {
00339         return 255 - i[d];
00340     }
00341 };
00342 
00343 template <class ImageType>
00344 HuginBase::CPVector _getVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,ImageType& image, vigra::BImage& mask, const unsigned int nrLines)
00345 {
00346     HuginBase::CPVector verticalLines;
00347     HuginBase::CPVector detectedLines;
00348     const HuginBase::SrcPanoImage& srcImage=pano.getImage(imgNr);
00349     bool needsRemap=srcImage.getProjection()!=HuginBase::SrcPanoImage::RECTILINEAR;
00350     double roll=(needsRemap?0:srcImage.getRoll());
00351     double size_factor=1.0;
00352     HuginBase::SrcPanoImage remappedImage;
00353     HuginBase::PanoramaOptions opts;
00354     vigra::BImage* edge;
00355     if(!needsRemap)
00356     {
00357         //rectilinear image can be used as is
00358         //detect edges
00359         edge=detectEdges(image,2,4,MAX_RESIZE_DIM,size_factor);
00360     }
00361     else
00362     {
00363         //remap all other image to equirectangular
00364         //create temporary SrcPanoImage, set appropriate image variables
00365         remappedImage=pano.getSrcImage(imgNr);
00366         remappedImage.setYaw(0);
00367         remappedImage.setPitch(0);
00368         remappedImage.setX(0);
00369         remappedImage.setY(0);
00370         remappedImage.setZ(0);
00371         remappedImage.setExposureValue(0);
00372         remappedImage.setEMoRParams(std::vector<float>(5, 0.0));
00373         remappedImage.deleteAllMasks();
00374         remappedImage.setActive(true);
00375         //create PanoramaOptions for remapping of image
00376         opts.setProjection(HuginBase::PanoramaOptions::EQUIRECTANGULAR);
00377         opts.setWidth(MAX_RESIZE_DIM);
00378         opts.outputExposureValue=0;
00379         //calculate output canvas size
00380         HuginBase::Panorama tempPano;
00381         tempPano.addImage(remappedImage);
00382         tempPano.setOptions(opts);
00383 
00384         HuginBase::CalculateFitPanorama fitPano(tempPano);
00385         fitPano.run();
00386         opts.setHFOV(fitPano.getResultHorizontalFOV());
00387         opts.setHeight(hugin_utils::roundi(fitPano.getResultHeight()));
00388         tempPano.setOptions(opts);
00389 
00390         //finally remap image
00391         HuginBase::Nona::RemappedPanoImage<ImageType,vigra::BImage>* remapped=new HuginBase::Nona::RemappedPanoImage<ImageType,vigra::BImage>;
00392         AppBase::ProgressDisplay* progress=new AppBase::DummyProgressDisplay();
00393         remapped->setPanoImage(remappedImage,opts,opts.getROI());
00394         if (mask.size().area() > 0)
00395         {
00396             remapped->remapImage(vigra::srcImageRange(image), vigra::srcImage(mask), vigra_ext::INTERP_CUBIC, progress);
00397         }
00398         else
00399         {
00400             remapped->remapImage(vigra::srcImageRange(image), vigra_ext::INTERP_CUBIC, progress);
00401         };
00402         ImageType remappedBitmap=remapped->m_image;
00403         mask = remapped->m_mask;
00404         //detect edges
00405         edge=detectEdges(remappedBitmap,2,4,std::max(remappedBitmap.width(),remappedBitmap.height())+10,size_factor);
00406         delete remapped;
00407         delete progress;
00408     };
00409     // ignore all edges outside of masked areas
00410     if (mask.size().area() > 0)
00411     {
00412         vigra::initImageIf(vigra::destImageRange(*edge), vigra::srcImage(mask, InvertedMaskAccessor()), vigra::UInt8(255));
00413     };
00414     //detect lines
00415     //we need the focal length
00416     double focalLength=srcImage.getExifFocalLength();
00417     if(focalLength==0)
00418     {
00419         focalLength=HuginBase::SrcPanoImage::calcFocalLength(
00420             srcImage.getProjection(),srcImage.getHFOV(),srcImage.getCropFactor(),srcImage.getSize());
00421     };
00422     Lines foundLines=findLines(*edge,0.05,focalLength,srcImage.getCropFactor());
00423     delete edge;
00424     //filter results
00425     VerticalLineVector filteredLines=FilterLines(foundLines,roll);
00426     //create control points
00427     if(!filteredLines.empty())
00428     {
00429         //we need to transform the coordinates to image coordinates because the detection
00430         //worked on smaller images or in remapped image
00431         HuginBase::PTools::Transform transform;
00432         if(needsRemap)
00433         {
00434             transform.createTransform(remappedImage,opts);
00435         };
00436         for(size_t i=0; i<filteredLines.size(); i++)
00437         {
00438             HuginBase::ControlPoint cp;
00439             cp.image1Nr=0;
00440             cp.image2Nr=0;
00441             cp.mode=HuginBase::ControlPoint::X;
00442             if(!needsRemap)
00443             {
00444                 cp.x1 = filteredLines[i].GetStart().x*size_factor;
00445                 cp.y1 = filteredLines[i].GetStart().y*size_factor;
00446                 cp.x2 = filteredLines[i].GetEnd().x*size_factor;
00447                 cp.y2 = filteredLines[i].GetEnd().y*size_factor;
00448             }
00449             else
00450             {
00451                 double xout;
00452                 double yout;
00453                 if(!transform.transformImgCoord(xout,yout,filteredLines[i].GetStart().x,filteredLines[i].GetStart().y))
00454                 {
00455                     continue;
00456                 };
00457                 cp.x1=xout;
00458                 cp.y1=yout;
00459                 if(!transform.transformImgCoord(xout,yout,filteredLines[i].GetEnd().x,filteredLines[i].GetEnd().y))
00460                 {
00461                     continue;
00462                 };
00463                 cp.x2=xout;
00464                 cp.y2=yout;
00465             };
00466             if(cp.x1>=0 && cp.x1<srcImage.getWidth() && cp.y1>=0 && cp.y1<srcImage.getHeight() &&
00467                cp.x2>=0 && cp.x2<srcImage.getWidth() && cp.y2>=0 && cp.y2<srcImage.getHeight())
00468             {
00469                 detectedLines.push_back(cp);
00470             };
00471         };
00472         //now a final check of the found vertical lines
00473         //we optimize the pano with a single image and disregard vertical lines with bigger errors
00474         //we need at least 2 lines
00475         if(detectedLines.size()>1)
00476         {
00477             HuginBase::Panorama tempPano;
00478             HuginBase::SrcPanoImage tempImage=pano.getSrcImage(imgNr);
00479             tempImage.setYaw(0);
00480             tempImage.setPitch(0);
00481             tempImage.setRoll(0);
00482             tempImage.setX(0);
00483             tempImage.setY(0);
00484             tempImage.setZ(0);
00485             tempPano.addImage(tempImage);
00486             for(size_t i=0; i<detectedLines.size(); i++)
00487             {
00488                 tempPano.addCtrlPoint(detectedLines[i]);
00489             };
00490             HuginBase::PanoramaOptions opt2;
00491             opt2.setProjection(HuginBase::PanoramaOptions::EQUIRECTANGULAR);
00492             tempPano.setOptions(opt2);
00493             HuginBase::OptimizeVector optVec;
00494             std::set<std::string> imgopt;
00495             imgopt.insert("p");
00496             imgopt.insert("r");
00497             optVec.push_back(imgopt);
00498             tempPano.setOptimizeVector(optVec);
00499             // ARGH the panotools optimizer uses global variables is not reentrant
00500 #pragma omp critical
00501             {
00502                 HuginBase::PTools::optimize(tempPano);
00503             }
00504             //first filter stage
00505             //we disregard all lines with big error
00506             //calculate statistic and determine limit
00507             double minError,maxError,mean,var;
00508             HuginBase::CalculateCPStatisticsError::calcCtrlPntsErrorStats(tempPano,minError,maxError,mean,var);
00509             detectedLines=tempPano.getCtrlPoints();
00510             double limit=mean+sqrt(var);
00511             maxError=0;
00512             for(int i=detectedLines.size()-1; i>=0; i--)
00513             {
00514                 if(detectedLines[i].error>limit)
00515                 {
00516                     detectedLines.erase(detectedLines.begin()+i);
00517                 }
00518                 else
00519                 {
00520                     //we need the max error of the remaining lines for the next step
00521                     maxError=std::max(detectedLines[i].error,maxError);
00522                 };
00523             };
00524             if(!detectedLines.empty() && maxError>0) //security check, should never be false
00525             {
00526                 //now keep only the best nrLines lines
00527                 //we are using error and line length as figure of merrit
00528                 for(size_t i=0;i<detectedLines.size();i++)
00529                 {
00530                     double length=sqrt(hugin_utils::sqr(detectedLines[i].x2-detectedLines[i].x1)+hugin_utils::sqr(detectedLines[i].y2-detectedLines[i].y1));
00531                     //calculate number of merrit
00532                     detectedLines[i].error=detectedLines[i].error/maxError+(1.0-std::min(length,500.0)/500.0);
00533                 };
00534                 std::sort(detectedLines.begin(),detectedLines.end(),SortByError);
00535                 //only save best nrLines control points
00536                 for(size_t i=0;i<detectedLines.size() && i<nrLines; i++)
00537                 {
00538                     HuginBase::ControlPoint cp=detectedLines[i];
00539                     cp.image1Nr=imgNr;
00540                     cp.image2Nr=imgNr;
00541                     cp.error=0;
00542                     verticalLines.push_back(cp);
00543                 };
00544             };
00545         }
00546         else
00547         {
00548             //if only one line was detected we do a special check
00549             //the allow deviation between line and roll angle is checked more narrow than in the first check
00550             if(detectedLines.size()==1)
00551             {
00552                 vigra::Diff2D diff((double)detectedLines[0].x2-detectedLines[0].x1,(double)detectedLines[0].y2-detectedLines[0].y1);
00553                 if(std::abs((diff.x*cos(DEG_TO_RAD(roll))+diff.y*sin(DEG_TO_RAD(roll)))/diff.magnitude())<0.05)
00554                 {
00555                     HuginBase::ControlPoint cp=detectedLines[0];
00556                     cp.image1Nr=imgNr;
00557                     cp.image2Nr=imgNr;
00558                     cp.error=0;
00559                     verticalLines.push_back(cp);
00560                 };
00561             };
00562         };
00563     };
00564     return verticalLines;
00565 };
00566 
00567 HuginBase::CPVector GetVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,vigra::UInt8RGBImage& image, vigra::BImage& mask, const unsigned int nrLines)
00568 {
00569     return _getVerticalLines(pano, imgNr, image, mask, nrLines);
00570 };
00571 
00572 HuginBase::CPVector GetVerticalLines(const HuginBase::Panorama& pano,const unsigned int imgNr,vigra::BImage& image, vigra::BImage& mask, const unsigned int nrLines)
00573 {
00574     return _getVerticalLines(pano, imgNr, image, mask, nrLines);
00575 };
00576 
00577 }; //namespace

Generated on 23 May 2018 for Hugintrunk by  doxygen 1.4.7