lsd.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                           License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
00014 // Third party copyrights are property of their respective owners.
00015 //
00016 // Redistribution and use in source and binary forms, with or without modification,
00017 // are permitted provided that the following conditions are met:
00018 //
00019 //   * Redistributions of source code must retain the above copyright notice,
00020 //     this list of conditions and the following disclaimer.
00021 //
00022 //   * Redistributions in binary form must reproduce the above copyright notice,
00023 //     this list of conditions and the following disclaimer in the documentation
00024 //     and/or other materials provided with the distribution.
00025 //
00026 //   * The name of the copyright holders may not be used to endorse or promote products
00027 //     derived from this software without specific prior written permission.
00028 //
00029 // This software is provided by the copyright holders and contributors "as is" and
00030 // any express or implied warranties, including, but not limited to, the implied
00031 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00032 // In no event shall the Intel Corporation or contributors be liable for any direct,
00033 // indirect, incidental, special, exemplary, or consequential damages
00034 // (including, but not limited to, procurement of substitute goods or services;
00035 // loss of use, data, or profits; or business interruption) however caused
00036 // and on any theory of liability, whether in contract, strict liability,
00037 // or tort (including negligence or otherwise) arising in any way out of
00038 // the use of this software, even if advised of the possibility of such damage.
00039 //
00040 //M*/
00041 
00042 #include "precomp.hpp"
00043 #include <vector>
00044 
00045 /////////////////////////////////////////////////////////////////////////////////////////
00046 // Default LSD parameters
00047 // SIGMA_SCALE 0.6    - Sigma for Gaussian filter is computed as sigma = sigma_scale/scale.
00048 // QUANT       2.0    - Bound to the quantization error on the gradient norm.
00049 // ANG_TH      22.5   - Gradient angle tolerance in degrees.
00050 // LOG_EPS     0.0    - Detection threshold: -log10(NFA) > log_eps
00051 // DENSITY_TH  0.7    - Minimal density of region points in rectangle.
00052 // N_BINS      1024   - Number of bins in pseudo-ordering of gradient modulus.
00053 
00054 #define M_3_2_PI    (3 * CV_PI) / 2   // 3/2 pi
00055 #define M_2__PI     (2 * CV_PI)         // 2 pi
00056 
00057 #ifndef M_LN10
00058 #define M_LN10      2.30258509299404568402
00059 #endif
00060 
00061 #define NOTDEF      double(-1024.0) // Label for pixels with undefined gradient.
00062 
00063 #define NOTUSED     0   // Label for pixels not used in yet.
00064 #define USED        1   // Label for pixels already used in detection.
00065 
00066 #define RELATIVE_ERROR_FACTOR 100.0
00067 
00068 const double DEG_TO_RADS = CV_PI / 180;
00069 
00070 #define log_gamma(x) ((x)>15.0?log_gamma_windschitl(x):log_gamma_lanczos(x))
00071 
00072 struct edge
00073 {
00074     cv::Point  p;
00075     bool taken;
00076 };
00077 
00078 /////////////////////////////////////////////////////////////////////////////////////////
00079 
00080 inline double distSq(const double x1, const double y1,
00081                      const double x2, const double y2)
00082 {
00083     return (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
00084 }
00085 
00086 inline double dist(const double x1, const double y1,
00087                    const double x2, const double y2)
00088 {
00089     return sqrt(distSq(x1, y1, x2, y2));
00090 }
00091 
00092 // Signed angle difference
00093 inline double angle_diff_signed(const double& a, const double& b)
00094 {
00095     double diff = a - b;
00096     while(diff <= -CV_PI) diff += M_2__PI;
00097     while(diff >   CV_PI) diff -= M_2__PI;
00098     return diff;
00099 }
00100 
00101 // Absolute value angle difference
00102 inline double angle_diff(const double& a, const double& b)
00103 {
00104     return std::fabs(angle_diff_signed(a, b));
00105 }
00106 
00107 // Compare doubles by relative error.
00108 inline bool double_equal(const double& a, const double& b)
00109 {
00110     // trivial case
00111     if(a == b) return true;
00112 
00113     double abs_diff = fabs(a - b);
00114     double aa = fabs(a);
00115     double bb = fabs(b);
00116     double abs_max = (aa > bb)? aa : bb;
00117 
00118     if(abs_max < DBL_MIN) abs_max = DBL_MIN;
00119 
00120     return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
00121 }
00122 
00123 inline bool AsmallerB_XoverY(const edge& a, const edge& b)
00124 {
00125     if (a.p.x == b.p.x) return a.p.y < b.p.y;
00126     else return a.p.x < b.p.x;
00127 }
00128 
00129 /**
00130  *   Computes the natural logarithm of the absolute value of
00131  *   the gamma function of x using Windschitl method.
00132  *   See http://www.rskey.org/gamma.htm
00133  */
00134 inline double log_gamma_windschitl(const double& x)
00135 {
00136     return 0.918938533204673 + (x-0.5)*log(x) - x
00137          + 0.5*x*log(x*sinh(1/x) + 1/(810.0*pow(x, 6.0)));
00138 }
00139 
00140 /**
00141  *   Computes the natural logarithm of the absolute value of
00142  *   the gamma function of x using the Lanczos approximation.
00143  *   See http://www.rskey.org/gamma.htm
00144  */
00145 inline double log_gamma_lanczos(const double& x)
00146 {
00147     static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
00148                          8687.24529705, 1168.92649479, 83.8676043424,
00149                          2.50662827511 };
00150     double a = (x + 0.5) * log(x + 5.5) - (x + 5.5);
00151     double b = 0;
00152     for(int n = 0; n < 7; ++n)
00153     {
00154         a -= log(x + double(n));
00155         b += q[n] * pow(x, double(n));
00156     }
00157     return a + log(b);
00158 }
00159 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
00160 
00161 namespace cv{
00162 
00163 class LineSegmentDetectorImpl : public LineSegmentDetector
00164 {
00165 public:
00166 
00167 /**
00168  * Create a LineSegmentDetectorImpl object. Specifying scale, number of subdivisions for the image, should the lines be refined and other constants as follows:
00169  *
00170  * @param _refine       How should the lines found be refined?
00171  *                      LSD_REFINE_NONE - No refinement applied.
00172  *                      LSD_REFINE_STD  - Standard refinement is applied. E.g. breaking arches into smaller line approximations.
00173  *                      LSD_REFINE_ADV  - Advanced refinement. Number of false alarms is calculated,
00174  *                                    lines are refined through increase of precision, decrement in size, etc.
00175  * @param _scale        The scale of the image that will be used to find the lines. Range (0..1].
00176  * @param _sigma_scale  Sigma for Gaussian filter is computed as sigma = _sigma_scale/_scale.
00177  * @param _quant        Bound to the quantization error on the gradient norm.
00178  * @param _ang_th       Gradient angle tolerance in degrees.
00179  * @param _log_eps      Detection threshold: -log10(NFA) > _log_eps
00180  * @param _density_th   Minimal density of aligned region points in rectangle.
00181  * @param _n_bins       Number of bins in pseudo-ordering of gradient modulus.
00182  */
00183     LineSegmentDetectorImpl(int _refine = LSD_REFINE_STD, double _scale = 0.8,
00184         double _sigma_scale = 0.6, double _quant = 2.0, double _ang_th = 22.5,
00185         double _log_eps = 0, double _density_th = 0.7, int _n_bins = 1024);
00186 
00187 /**
00188  * Detect lines in the input image.
00189  *
00190  * @param _image    A grayscale(CV_8UC1) input image.
00191  *                  If only a roi needs to be selected, use
00192  *                  lsd_ptr->detect(image(roi), ..., lines);
00193  *                  lines += Scalar(roi.x, roi.y, roi.x, roi.y);
00194  * @param _lines    Return: A vector of Vec4i or Vec4f elements specifying the beginning and ending point of a line.
00195  *                          Where Vec4i/Vec4f is (x1, y1, x2, y2), point 1 is the start, point 2 - end.
00196  *                          Returned lines are strictly oriented depending on the gradient.
00197  * @param width     Return: Vector of widths of the regions, where the lines are found. E.g. Width of line.
00198  * @param prec      Return: Vector of precisions with which the lines are found.
00199  * @param nfa       Return: Vector containing number of false alarms in the line region, with precision of 10%.
00200  *                          The bigger the value, logarithmically better the detection.
00201  *                              * -1 corresponds to 10 mean false alarms
00202  *                              * 0 corresponds to 1 mean false alarm
00203  *                              * 1 corresponds to 0.1 mean false alarms
00204  *                          This vector will be calculated _only_ when the objects type is REFINE_ADV
00205  */
00206     void detect(InputArray _image, OutputArray _lines,
00207                 OutputArray width = noArray(), OutputArray prec = noArray(),
00208                 OutputArray nfa = noArray());
00209 
00210 /**
00211  * Draw lines on the given canvas.
00212  *
00213  * @param image     The image, where lines will be drawn.
00214  *                  Should have the size of the image, where the lines were found
00215  * @param lines     The lines that need to be drawn
00216  */
00217     void drawSegments(InputOutputArray _image, InputArray lines);
00218 
00219 /**
00220  * Draw both vectors on the image canvas. Uses blue for lines 1 and red for lines 2.
00221  *
00222  * @param size      The size of the image, where lines1 and lines2 were found.
00223  * @param lines1    The first lines that need to be drawn. Color - Blue.
00224  * @param lines2    The second lines that need to be drawn. Color - Red.
00225  * @param image     An optional image, where lines will be drawn.
00226  *                  Should have the size of the image, where the lines were found
00227  * @return          The number of mismatching pixels between lines1 and lines2.
00228  */
00229     int compareSegments(const Size& size, InputArray lines1, InputArray lines2, InputOutputArray _image = noArray());
00230 
00231 private:
00232     Mat image;
00233     Mat_<double> scaled_image;
00234     double *scaled_image_data;
00235     Mat_<double> angles;     // in rads
00236     double *angles_data;
00237     Mat_<double> modgrad;
00238     double *modgrad_data;
00239     Mat_<uchar> used;
00240 
00241     int img_width;
00242     int img_height;
00243     double LOG_NT;
00244 
00245     bool w_needed;
00246     bool p_needed;
00247     bool n_needed;
00248 
00249     const double SCALE;
00250     const int doRefine;
00251     const double SIGMA_SCALE;
00252     const double QUANT;
00253     const double ANG_TH;
00254     const double LOG_EPS;
00255     const double DENSITY_TH;
00256     const int N_BINS;
00257 
00258     struct RegionPoint {
00259         int x;
00260         int y;
00261         uchar* used;
00262         double angle;
00263         double modgrad;
00264     };
00265 
00266 
00267     struct coorlist
00268     {
00269         Point2i p;
00270         struct coorlist* next;
00271     };
00272 
00273     struct rect
00274     {
00275         double x1, y1, x2, y2;    // first and second point of the line segment
00276         double width;             // rectangle width
00277         double x, y;              // center of the rectangle
00278         double theta;             // angle
00279         double dx,dy;             // (dx,dy) is vector oriented as the line segment
00280         double prec;              // tolerance angle
00281         double p;                 // probability of a point with angle within 'prec'
00282     };
00283 
00284     LineSegmentDetectorImpl& operator= (const LineSegmentDetectorImpl&); // to quiet MSVC
00285 
00286 /**
00287  * Detect lines in the whole input image.
00288  *
00289  * @param lines         Return: A vector of Vec4f elements specifying the beginning and ending point of a line.
00290  *                              Where Vec4f is (x1, y1, x2, y2), point 1 is the start, point 2 - end.
00291  *                              Returned lines are strictly oriented depending on the gradient.
00292  * @param widths        Return: Vector of widths of the regions, where the lines are found. E.g. Width of line.
00293  * @param precisions    Return: Vector of precisions with which the lines are found.
00294  * @param nfas          Return: Vector containing number of false alarms in the line region, with precision of 10%.
00295  *                              The bigger the value, logarithmically better the detection.
00296  *                                  * -1 corresponds to 10 mean false alarms
00297  *                                  * 0 corresponds to 1 mean false alarm
00298  *                                  * 1 corresponds to 0.1 mean false alarms
00299  */
00300     void flsd(std::vector<Vec4f>& lines,
00301               std::vector<double>& widths, std::vector<double>& precisions,
00302               std::vector<double>& nfas);
00303 
00304 /**
00305  * Finds the angles and the gradients of the image. Generates a list of pseudo ordered points.
00306  *
00307  * @param threshold The minimum value of the angle that is considered defined, otherwise NOTDEF
00308  * @param n_bins    The number of bins with which gradients are ordered by, using bucket sort.
00309  * @param list      Return: Vector of coordinate points that are pseudo ordered by magnitude.
00310  *                  Pixels would be ordered by norm value, up to a precision given by max_grad/n_bins.
00311  */
00312     void ll_angle(const double& threshold, const unsigned int& n_bins, std::vector<coorlist>& list);
00313 
00314 /**
00315  * Grow a region starting from point s with a defined precision,
00316  * returning the containing points size and the angle of the gradients.
00317  *
00318  * @param s         Starting point for the region.
00319  * @param reg       Return: Vector of points, that are part of the region
00320  * @param reg_size  Return: The size of the region.
00321  * @param reg_angle Return: The mean angle of the region.
00322  * @param prec      The precision by which each region angle should be aligned to the mean.
00323  */
00324     void region_grow(const Point2i& s, std::vector<RegionPoint>& reg,
00325                      int& reg_size, double& reg_angle, const double& prec);
00326 
00327 /**
00328  * Finds the bounding rotated rectangle of a region.
00329  *
00330  * @param reg       The region of points, from which the rectangle to be constructed from.
00331  * @param reg_size  The number of points in the region.
00332  * @param reg_angle The mean angle of the region.
00333  * @param prec      The precision by which points were found.
00334  * @param p         Probability of a point with angle within 'prec'.
00335  * @param rec       Return: The generated rectangle.
00336  */
00337     void region2rect(const std::vector<RegionPoint>& reg, const int reg_size, const double reg_angle,
00338                      const double prec, const double p, rect& rec) const;
00339 
00340 /**
00341  * Compute region's angle as the principal inertia axis of the region.
00342  * @return          Regions angle.
00343  */
00344     double get_theta(const std::vector<RegionPoint>& reg, const int& reg_size, const double& x,
00345                      const double& y, const double& reg_angle, const double& prec) const;
00346 
00347 /**
00348  * An estimation of the angle tolerance is performed by the standard deviation of the angle at points
00349  * near the region's starting point. Then, a new region is grown starting from the same point, but using the
00350  * estimated angle tolerance. If this fails to produce a rectangle with the right density of region points,
00351  * 'reduce_region_radius' is called to try to satisfy this condition.
00352  */
00353     bool refine(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
00354                 const double prec, double p, rect& rec, const double& density_th);
00355 
00356 /**
00357  * Reduce the region size, by elimination the points far from the starting point, until that leads to
00358  * rectangle with the right density of region points or to discard the region if too small.
00359  */
00360     bool reduce_region_radius(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
00361                 const double prec, double p, rect& rec, double density, const double& density_th);
00362 
00363 /**
00364  * Try some rectangles variations to improve NFA value. Only if the rectangle is not meaningful (i.e., log_nfa <= log_eps).
00365  * @return      The new NFA value.
00366  */
00367     double rect_improve(rect& rec) const;
00368 
00369 /**
00370  * Calculates the number of correctly aligned points within the rectangle.
00371  * @return      The new NFA value.
00372  */
00373     double rect_nfa(const rect& rec) const;
00374 
00375 /**
00376  * Computes the NFA values based on the total number of points, points that agree.
00377  * n, k, p are the binomial parameters.
00378  * @return      The new NFA value.
00379  */
00380     double nfa(const int& n, const int& k, const double& p) const;
00381 
00382 /**
00383  * Is the point at place 'address' aligned to angle theta, up to precision 'prec'?
00384  * @return      Whether the point is aligned.
00385  */
00386     bool isAligned(const int& address, const double& theta, const double& prec) const;
00387 };
00388 
00389 /////////////////////////////////////////////////////////////////////////////////////////
00390 
00391 CV_EXPORTS Ptr<LineSegmentDetector> createLineSegmentDetector(
00392         int _refine, double _scale, double _sigma_scale, double _quant, double _ang_th,
00393         double _log_eps, double _density_th, int _n_bins)
00394 {
00395     return makePtr<LineSegmentDetectorImpl>(
00396             _refine, _scale, _sigma_scale, _quant, _ang_th,
00397             _log_eps, _density_th, _n_bins);
00398 }
00399 
00400 /////////////////////////////////////////////////////////////////////////////////////////
00401 
00402 LineSegmentDetectorImpl::LineSegmentDetectorImpl(int _refine, double _scale, double _sigma_scale, double _quant,
00403         double _ang_th, double _log_eps, double _density_th, int _n_bins)
00404         :SCALE(_scale), doRefine(_refine), SIGMA_SCALE(_sigma_scale), QUANT(_quant),
00405         ANG_TH(_ang_th), LOG_EPS(_log_eps), DENSITY_TH(_density_th), N_BINS(_n_bins)
00406 {
00407     CV_Assert(_scale > 0 && _sigma_scale > 0 && _quant >= 0 &&
00408               _ang_th > 0 && _ang_th < 180 && _density_th >= 0 && _density_th < 1 &&
00409               _n_bins > 0);
00410 }
00411 
00412 void LineSegmentDetectorImpl::detect(InputArray _image, OutputArray _lines,
00413                 OutputArray _width, OutputArray _prec, OutputArray _nfa)
00414 {
00415     Mat_<double> img = _image.getMat();
00416     CV_Assert(!img.empty() && img.channels() == 1);
00417 
00418     // Convert image to double
00419     img.convertTo(image, CV_64FC1);
00420 
00421     std::vector<Vec4f> lines;
00422     std::vector<double> w, p, n;
00423     w_needed = _width.needed();
00424     p_needed = _prec.needed();
00425     if (doRefine < LSD_REFINE_ADV)
00426         n_needed = false;
00427     else
00428         n_needed = _nfa.needed();
00429 
00430     flsd(lines, w, p, n);
00431 
00432     Mat(lines).copyTo(_lines);
00433     if(w_needed) Mat(w).copyTo(_width);
00434     if(p_needed) Mat(p).copyTo(_prec);
00435     if(n_needed) Mat(n).copyTo(_nfa);
00436 }
00437 
00438 void LineSegmentDetectorImpl::flsd(std::vector<Vec4f>& lines,
00439     std::vector<double>& widths, std::vector<double>& precisions,
00440     std::vector<double>& nfas)
00441 {
00442     // Angle tolerance
00443     const double prec = CV_PI * ANG_TH / 180;
00444     const double p = ANG_TH / 180;
00445     const double rho = QUANT / sin(prec);    // gradient magnitude threshold
00446 
00447     std::vector<coorlist> list;
00448     if(SCALE != 1)
00449     {
00450         Mat gaussian_img;
00451         const double sigma = (SCALE < 1)?(SIGMA_SCALE / SCALE):(SIGMA_SCALE);
00452         const double sprec = 3;
00453         const unsigned int h =  (unsigned int)(ceil(sigma * sqrt(2 * sprec * log(10.0))));
00454         Size ksize(1 + 2 * h, 1 + 2 * h); // kernel size
00455         GaussianBlur(image, gaussian_img, ksize, sigma);
00456         // Scale image to needed size
00457         resize(gaussian_img, scaled_image, Size(), SCALE, SCALE);
00458         ll_angle(rho, N_BINS, list);
00459     }
00460     else
00461     {
00462         scaled_image = image;
00463         ll_angle(rho, N_BINS, list);
00464     }
00465 
00466     LOG_NT = 5 * (log10(double(img_width)) + log10(double(img_height))) / 2 + log10(11.0);
00467     const int min_reg_size = int(-LOG_NT/log10(p)); // minimal number of points in region that can give a meaningful event
00468 
00469     // // Initialize region only when needed
00470     // Mat region = Mat::zeros(scaled_image.size(), CV_8UC1);
00471     used = Mat_<uchar>::zeros(scaled_image.size()); // zeros = NOTUSED
00472     std::vector<RegionPoint> reg(img_width * img_height);
00473 
00474     // Search for line segments
00475     unsigned int ls_count = 0;
00476     for(size_t i = 0, list_size = list.size(); i < list_size; ++i)
00477     {
00478         unsigned int adx = list[i].p.x + list[i].p.y * img_width;
00479         if((used.ptr()[adx] == NOTUSED) && (angles_data[adx] != NOTDEF))
00480         {
00481             int reg_size;
00482             double reg_angle;
00483             region_grow(list[i].p, reg, reg_size, reg_angle, prec);
00484 
00485             // Ignore small regions
00486             if(reg_size < min_reg_size) { continue; }
00487 
00488             // Construct rectangular approximation for the region
00489             rect rec;
00490             region2rect(reg, reg_size, reg_angle, prec, p, rec);
00491 
00492             double log_nfa = -1;
00493             if(doRefine > LSD_REFINE_NONE)
00494             {
00495                 // At least REFINE_STANDARD lvl.
00496                 if(!refine(reg, reg_size, reg_angle, prec, p, rec, DENSITY_TH)) { continue; }
00497 
00498                 if(doRefine >= LSD_REFINE_ADV)
00499                 {
00500                     // Compute NFA
00501                     log_nfa = rect_improve(rec);
00502                     if(log_nfa <= LOG_EPS) { continue; }
00503                 }
00504             }
00505             // Found new line
00506             ++ls_count;
00507 
00508             // Add the offset
00509             rec.x1 += 0.5; rec.y1 += 0.5;
00510             rec.x2 += 0.5; rec.y2 += 0.5;
00511 
00512             // scale the result values if a sub-sampling was performed
00513             if(SCALE != 1)
00514             {
00515                 rec.x1 /= SCALE; rec.y1 /= SCALE;
00516                 rec.x2 /= SCALE; rec.y2 /= SCALE;
00517                 rec.width /= SCALE;
00518             }
00519 
00520             //Store the relevant data
00521             lines.push_back(Vec4f(float(rec.x1), float(rec.y1), float(rec.x2), float(rec.y2)));
00522             if(w_needed) widths.push_back(rec.width);
00523             if(p_needed) precisions.push_back(rec.p);
00524             if(n_needed && doRefine >= LSD_REFINE_ADV) nfas.push_back(log_nfa);
00525 
00526 
00527             // //Add the linesID to the region on the image
00528             // for(unsigned int el = 0; el < reg_size; el++)
00529             // {
00530             //     region.data[reg[i].x + reg[i].y * width] = ls_count;
00531             // }
00532         }
00533     }
00534 }
00535 
00536 void LineSegmentDetectorImpl::ll_angle(const double& threshold,
00537                                    const unsigned int& n_bins,
00538                                    std::vector<coorlist>& list)
00539 {
00540     //Initialize data
00541     angles = Mat_<double>(scaled_image.size());
00542     modgrad = Mat_<double>(scaled_image.size());
00543 
00544     angles_data = angles.ptr<double>(0);
00545     modgrad_data = modgrad.ptr<double>(0);
00546     scaled_image_data = scaled_image.ptr<double>(0);
00547 
00548     img_width = scaled_image.cols;
00549     img_height = scaled_image.rows;
00550 
00551     // Undefined the down and right boundaries
00552     angles.row(img_height - 1).setTo(NOTDEF);
00553     angles.col(img_width - 1).setTo(NOTDEF);
00554 
00555     // Computing gradient for remaining pixels
00556     CV_Assert(scaled_image.isContinuous() &&
00557               modgrad.isContinuous() &&
00558               angles.isContinuous());   // Accessing image data linearly
00559 
00560     double max_grad = -1;
00561     for(int y = 0; y < img_height - 1; ++y)
00562     {
00563         for(int addr = y * img_width, addr_end = addr + img_width - 1; addr < addr_end; ++addr)
00564         {
00565             double DA = scaled_image_data[addr + img_width + 1] - scaled_image_data[addr];
00566             double BC = scaled_image_data[addr + 1] - scaled_image_data[addr + img_width];
00567             double gx = DA + BC;    // gradient x component
00568             double gy = DA - BC;    // gradient y component
00569             double norm = std::sqrt((gx * gx + gy * gy) / 4); // gradient norm
00570 
00571             modgrad_data[addr] = norm;    // store gradient
00572 
00573             if (norm <= threshold)  // norm too small, gradient no defined
00574             {
00575                 angles_data[addr] = NOTDEF;
00576             }
00577             else
00578             {
00579                 angles_data[addr] = fastAtan2(float(gx), float(-gy)) * DEG_TO_RADS;  // gradient angle computation
00580                 if (norm > max_grad) { max_grad = norm; }
00581             }
00582 
00583         }
00584     }
00585 
00586     // Compute histogram of gradient values
00587     list = std::vector<coorlist>(img_width * img_height);
00588     std::vector<coorlist*> range_s(n_bins);
00589     std::vector<coorlist*> range_e(n_bins);
00590     unsigned int count = 0;
00591     double bin_coef = (max_grad > 0) ? double(n_bins - 1) / max_grad : 0; // If all image is smooth, max_grad <= 0
00592 
00593     for(int y = 0; y < img_height - 1; ++y)
00594     {
00595         const double* norm = modgrad_data + y * img_width;
00596         for(int x = 0; x < img_width - 1; ++x, ++norm)
00597         {
00598             // Store the point in the right bin according to its norm
00599             int i = int((*norm) * bin_coef);
00600             if(!range_e[i])
00601             {
00602                 range_e[i] = range_s[i] = &list[count];
00603                 ++count;
00604             }
00605             else
00606             {
00607                 range_e[i]->next = &list[count];
00608                 range_e[i] = &list[count];
00609                 ++count;
00610             }
00611             range_e[i]->p = Point(x, y);
00612             range_e[i]->next = 0;
00613         }
00614     }
00615 
00616     // Sort
00617     int idx = n_bins - 1;
00618     for(;idx > 0 && !range_s[idx]; --idx);
00619     coorlist* start = range_s[idx];
00620     coorlist* end = range_e[idx];
00621     if(start)
00622     {
00623         while(idx > 0)
00624         {
00625             --idx;
00626             if(range_s[idx])
00627             {
00628                 end->next = range_s[idx];
00629                 end = range_e[idx];
00630             }
00631         }
00632     }
00633 }
00634 
00635 void LineSegmentDetectorImpl::region_grow(const Point2i& s, std::vector<RegionPoint>& reg,
00636                                       int& reg_size, double& reg_angle, const double& prec)
00637 {
00638     // Point to this region
00639     reg_size = 1;
00640     reg[0].x = s.x;
00641     reg[0].y = s.y;
00642     int addr = s.x + s.y * img_width;
00643     reg[0].used = used.ptr() + addr;
00644     reg_angle = angles_data[addr];
00645     reg[0].angle = reg_angle;
00646     reg[0].modgrad = modgrad_data[addr];
00647 
00648     float sumdx = float(std::cos(reg_angle));
00649     float sumdy = float(std::sin(reg_angle));
00650     *reg[0].used = USED;
00651 
00652     //Try neighboring regions
00653     for(int i = 0; i < reg_size; ++i)
00654     {
00655         const RegionPoint& rpoint = reg[i];
00656         int xx_min = std::max(rpoint.x - 1, 0), xx_max = std::min(rpoint.x + 1, img_width - 1);
00657         int yy_min = std::max(rpoint.y - 1, 0), yy_max = std::min(rpoint.y + 1, img_height - 1);
00658         for(int yy = yy_min; yy <= yy_max; ++yy)
00659         {
00660             int c_addr = xx_min + yy * img_width;
00661             for(int xx = xx_min; xx <= xx_max; ++xx, ++c_addr)
00662             {
00663                 if((used.ptr()[c_addr] != USED) &&
00664                    (isAligned(c_addr, reg_angle, prec)))
00665                 {
00666                     // Add point
00667                     used.ptr()[c_addr] = USED;
00668                     RegionPoint& region_point = reg[reg_size];
00669                     region_point.x = xx;
00670                     region_point.y = yy;
00671                     region_point.used = &(used.ptr()[c_addr]);
00672                     region_point.modgrad = modgrad_data[c_addr];
00673                     const double& angle = angles_data[c_addr];
00674                     region_point.angle = angle;
00675                     ++reg_size;
00676 
00677                     // Update region's angle
00678                     sumdx += cos(float(angle));
00679                     sumdy += sin(float(angle));
00680                     // reg_angle is used in the isAligned, so it needs to be updates?
00681                     reg_angle = fastAtan2(sumdy, sumdx) * DEG_TO_RADS;
00682                 }
00683             }
00684         }
00685     }
00686 }
00687 
00688 void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg, const int reg_size,
00689                                       const double reg_angle, const double prec, const double p, rect& rec) const
00690 {
00691     double x = 0, y = 0, sum = 0;
00692     for(int i = 0; i < reg_size; ++i)
00693     {
00694         const RegionPoint& pnt = reg[i];
00695         const double& weight = pnt.modgrad;
00696         x += double(pnt.x) * weight;
00697         y += double(pnt.y) * weight;
00698         sum += weight;
00699     }
00700 
00701     // Weighted sum must differ from 0
00702     CV_Assert(sum > 0);
00703 
00704     x /= sum;
00705     y /= sum;
00706 
00707     double theta = get_theta(reg, reg_size, x, y, reg_angle, prec);
00708 
00709     // Find length and width
00710     double dx = cos(theta);
00711     double dy = sin(theta);
00712     double l_min = 0, l_max = 0, w_min = 0, w_max = 0;
00713 
00714     for(int i = 0; i < reg_size; ++i)
00715     {
00716         double regdx = double(reg[i].x) - x;
00717         double regdy = double(reg[i].y) - y;
00718 
00719         double l = regdx * dx + regdy * dy;
00720         double w = -regdx * dy + regdy * dx;
00721 
00722         if(l > l_max) l_max = l;
00723         else if(l < l_min) l_min = l;
00724         if(w > w_max) w_max = w;
00725         else if(w < w_min) w_min = w;
00726     }
00727 
00728     // Store values
00729     rec.x1 = x + l_min * dx;
00730     rec.y1 = y + l_min * dy;
00731     rec.x2 = x + l_max * dx;
00732     rec.y2 = y + l_max * dy;
00733     rec.width = w_max - w_min;
00734     rec.x = x;
00735     rec.y = y;
00736     rec.theta = theta;
00737     rec.dx = dx;
00738     rec.dy = dy;
00739     rec.prec = prec;
00740     rec.p = p;
00741 
00742     // Min width of 1 pixel
00743     if(rec.width < 1.0) rec.width = 1.0;
00744 }
00745 
00746 double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, const int& reg_size, const double& x,
00747                                       const double& y, const double& reg_angle, const double& prec) const
00748 {
00749     double Ixx = 0.0;
00750     double Iyy = 0.0;
00751     double Ixy = 0.0;
00752 
00753     // Compute inertia matrix
00754     for(int i = 0; i < reg_size; ++i)
00755     {
00756         const double& regx = reg[i].x;
00757         const double& regy = reg[i].y;
00758         const double& weight = reg[i].modgrad;
00759         double dx = regx - x;
00760         double dy = regy - y;
00761         Ixx += dy * dy * weight;
00762         Iyy += dx * dx * weight;
00763         Ixy -= dx * dy * weight;
00764     }
00765 
00766     // Check if inertia matrix is null
00767     CV_Assert(!(double_equal(Ixx, 0) && double_equal(Iyy, 0) && double_equal(Ixy, 0)));
00768 
00769     // Compute smallest eigenvalue
00770     double lambda = 0.5 * (Ixx + Iyy - sqrt((Ixx - Iyy) * (Ixx - Iyy) + 4.0 * Ixy * Ixy));
00771 
00772     // Compute angle
00773     double theta = (fabs(Ixx)>fabs(Iyy))?
00774                     double(fastAtan2(float(lambda - Ixx), float(Ixy))):
00775                     double(fastAtan2(float(Ixy), float(lambda - Iyy))); // in degs
00776     theta *= DEG_TO_RADS;
00777 
00778     // Correct angle by 180 deg if necessary
00779     if(angle_diff(theta, reg_angle) > prec) { theta += CV_PI; }
00780 
00781     return theta;
00782 }
00783 
00784 bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
00785                                  const double prec, double p, rect& rec, const double& density_th)
00786 {
00787     double density = double(reg_size) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
00788 
00789     if (density >= density_th) { return true; }
00790 
00791     // Try to reduce angle tolerance
00792     double xc = double(reg[0].x);
00793     double yc = double(reg[0].y);
00794     const double& ang_c = reg[0].angle;
00795     double sum = 0, s_sum = 0;
00796     int n = 0;
00797 
00798     for (int i = 0; i < reg_size; ++i)
00799     {
00800         *(reg[i].used) = NOTUSED;
00801         if (dist(xc, yc, reg[i].x, reg[i].y) < rec.width)
00802         {
00803             const double& angle = reg[i].angle;
00804             double ang_d = angle_diff_signed(angle, ang_c);
00805             sum += ang_d;
00806             s_sum += ang_d * ang_d;
00807             ++n;
00808         }
00809     }
00810     double mean_angle = sum / double(n);
00811     // 2 * standard deviation
00812     double tau = 2.0 * sqrt((s_sum - 2.0 * mean_angle * sum) / double(n) + mean_angle * mean_angle);
00813 
00814     // Try new region
00815     region_grow(Point(reg[0].x, reg[0].y), reg, reg_size, reg_angle, tau);
00816 
00817     if (reg_size < 2) { return false; }
00818 
00819     region2rect(reg, reg_size, reg_angle, prec, p, rec);
00820     density = double(reg_size) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
00821 
00822     if (density < density_th)
00823     {
00824         return reduce_region_radius(reg, reg_size, reg_angle, prec, p, rec, density, density_th);
00825     }
00826     else
00827     {
00828         return true;
00829     }
00830 }
00831 
00832 bool LineSegmentDetectorImpl::reduce_region_radius(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
00833                 const double prec, double p, rect& rec, double density, const double& density_th)
00834 {
00835     // Compute region's radius
00836     double xc = double(reg[0].x);
00837     double yc = double(reg[0].y);
00838     double radSq1 = distSq(xc, yc, rec.x1, rec.y1);
00839     double radSq2 = distSq(xc, yc, rec.x2, rec.y2);
00840     double radSq = radSq1 > radSq2 ? radSq1 : radSq2;
00841 
00842     while(density < density_th)
00843     {
00844         radSq *= 0.75*0.75; // Reduce region's radius to 75% of its value
00845         // Remove points from the region and update 'used' map
00846         for(int i = 0; i < reg_size; ++i)
00847         {
00848             if(distSq(xc, yc, double(reg[i].x), double(reg[i].y)) > radSq)
00849             {
00850                 // Remove point from the region
00851                 *(reg[i].used) = NOTUSED;
00852                 std::swap(reg[i], reg[reg_size - 1]);
00853                 --reg_size;
00854                 --i; // To avoid skipping one point
00855             }
00856         }
00857 
00858         if(reg_size < 2) { return false; }
00859 
00860         // Re-compute rectangle
00861         region2rect(reg, reg_size ,reg_angle, prec, p, rec);
00862 
00863         // Re-compute region points density
00864         density = double(reg_size) /
00865                   (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
00866     }
00867 
00868     return true;
00869 }
00870 
00871 double LineSegmentDetectorImpl::rect_improve(rect& rec) const
00872 {
00873     double delta = 0.5;
00874     double delta_2 = delta / 2.0;
00875 
00876     double log_nfa = rect_nfa(rec);
00877 
00878     if(log_nfa > LOG_EPS) return log_nfa; // Good rectangle
00879 
00880     // Try to improve
00881     // Finer precision
00882     rect r = rect(rec); // Copy
00883     for(int n = 0; n < 5; ++n)
00884     {
00885         r.p /= 2;
00886         r.prec = r.p * CV_PI;
00887         double log_nfa_new = rect_nfa(r);
00888         if(log_nfa_new > log_nfa)
00889         {
00890             log_nfa = log_nfa_new;
00891             rec = rect(r);
00892         }
00893     }
00894     if(log_nfa > LOG_EPS) return log_nfa;
00895 
00896     // Try to reduce width
00897     r = rect(rec);
00898     for(unsigned int n = 0; n < 5; ++n)
00899     {
00900         if((r.width - delta) >= 0.5)
00901         {
00902             r.width -= delta;
00903             double log_nfa_new = rect_nfa(r);
00904             if(log_nfa_new > log_nfa)
00905             {
00906                 rec = rect(r);
00907                 log_nfa = log_nfa_new;
00908             }
00909         }
00910     }
00911     if(log_nfa > LOG_EPS) return log_nfa;
00912 
00913     // Try to reduce one side of rectangle
00914     r = rect(rec);
00915     for(unsigned int n = 0; n < 5; ++n)
00916     {
00917         if((r.width - delta) >= 0.5)
00918         {
00919             r.x1 += -r.dy * delta_2;
00920             r.y1 +=  r.dx * delta_2;
00921             r.x2 += -r.dy * delta_2;
00922             r.y2 +=  r.dx * delta_2;
00923             r.width -= delta;
00924             double log_nfa_new = rect_nfa(r);
00925             if(log_nfa_new > log_nfa)
00926             {
00927                 rec = rect(r);
00928                 log_nfa = log_nfa_new;
00929             }
00930         }
00931     }
00932     if(log_nfa > LOG_EPS) return log_nfa;
00933 
00934     // Try to reduce other side of rectangle
00935     r = rect(rec);
00936     for(unsigned int n = 0; n < 5; ++n)
00937     {
00938         if((r.width - delta) >= 0.5)
00939         {
00940             r.x1 -= -r.dy * delta_2;
00941             r.y1 -=  r.dx * delta_2;
00942             r.x2 -= -r.dy * delta_2;
00943             r.y2 -=  r.dx * delta_2;
00944             r.width -= delta;
00945             double log_nfa_new = rect_nfa(r);
00946             if(log_nfa_new > log_nfa)
00947             {
00948                 rec = rect(r);
00949                 log_nfa = log_nfa_new;
00950             }
00951         }
00952     }
00953     if(log_nfa > LOG_EPS) return log_nfa;
00954 
00955     // Try finer precision
00956     r = rect(rec);
00957     for(unsigned int n = 0; n < 5; ++n)
00958     {
00959         if((r.width - delta) >= 0.5)
00960         {
00961             r.p /= 2;
00962             r.prec = r.p * CV_PI;
00963             double log_nfa_new = rect_nfa(r);
00964             if(log_nfa_new > log_nfa)
00965             {
00966                 rec = rect(r);
00967                 log_nfa = log_nfa_new;
00968             }
00969         }
00970     }
00971 
00972     return log_nfa;
00973 }
00974 
00975 double LineSegmentDetectorImpl::rect_nfa(const rect& rec) const
00976 {
00977     int total_pts = 0, alg_pts = 0;
00978     double half_width = rec.width / 2.0;
00979     double dyhw = rec.dy * half_width;
00980     double dxhw = rec.dx * half_width;
00981 
00982     std::vector<edge> ordered_x(4);
00983     edge* min_y = &ordered_x[0];
00984     edge* max_y = &ordered_x[0]; // Will be used for loop range
00985 
00986     ordered_x[0].p.x = int(rec.x1 - dyhw); ordered_x[0].p.y = int(rec.y1 + dxhw); ordered_x[0].taken = false;
00987     ordered_x[1].p.x = int(rec.x2 - dyhw); ordered_x[1].p.y = int(rec.y2 + dxhw); ordered_x[1].taken = false;
00988     ordered_x[2].p.x = int(rec.x2 + dyhw); ordered_x[2].p.y = int(rec.y2 - dxhw); ordered_x[2].taken = false;
00989     ordered_x[3].p.x = int(rec.x1 + dyhw); ordered_x[3].p.y = int(rec.y1 - dxhw); ordered_x[3].taken = false;
00990 
00991     std::sort(ordered_x.begin(), ordered_x.end(), AsmallerB_XoverY);
00992 
00993     // Find min y. And mark as taken. find max y.
00994     for(unsigned int i = 1; i < 4; ++i)
00995     {
00996         if(min_y->p.y > ordered_x[i].p.y) {min_y = &ordered_x[i]; }
00997         if(max_y->p.y < ordered_x[i].p.y) {max_y = &ordered_x[i]; }
00998     }
00999     min_y->taken = true;
01000 
01001     // Find leftmost untaken point;
01002     edge* leftmost = 0;
01003     for(unsigned int i = 0; i < 4; ++i)
01004     {
01005         if(!ordered_x[i].taken)
01006         {
01007             if(!leftmost) // if uninitialized
01008             {
01009                 leftmost = &ordered_x[i];
01010             }
01011             else if (leftmost->p.x > ordered_x[i].p.x)
01012             {
01013                 leftmost = &ordered_x[i];
01014             }
01015         }
01016     }
01017     leftmost->taken = true;
01018 
01019     // Find rightmost untaken point;
01020     edge* rightmost = 0;
01021     for(unsigned int i = 0; i < 4; ++i)
01022     {
01023         if(!ordered_x[i].taken)
01024         {
01025             if(!rightmost) // if uninitialized
01026             {
01027                 rightmost = &ordered_x[i];
01028             }
01029             else if (rightmost->p.x < ordered_x[i].p.x)
01030             {
01031                 rightmost = &ordered_x[i];
01032             }
01033         }
01034     }
01035     rightmost->taken = true;
01036 
01037     // Find last untaken point;
01038     edge* tailp = 0;
01039     for(unsigned int i = 0; i < 4; ++i)
01040     {
01041         if(!ordered_x[i].taken)
01042         {
01043             if(!tailp) // if uninitialized
01044             {
01045                 tailp = &ordered_x[i];
01046             }
01047             else if (tailp->p.x > ordered_x[i].p.x)
01048             {
01049                 tailp = &ordered_x[i];
01050             }
01051         }
01052     }
01053     tailp->taken = true;
01054 
01055     double flstep = (min_y->p.y != leftmost->p.y) ?
01056                     (min_y->p.x - leftmost->p.x) / (min_y->p.y - leftmost->p.y) : 0; //first left step
01057     double slstep = (leftmost->p.y != tailp->p.x) ?
01058                     (leftmost->p.x - tailp->p.x) / (leftmost->p.y - tailp->p.x) : 0; //second left step
01059 
01060     double frstep = (min_y->p.y != rightmost->p.y) ?
01061                     (min_y->p.x - rightmost->p.x) / (min_y->p.y - rightmost->p.y) : 0; //first right step
01062     double srstep = (rightmost->p.y != tailp->p.x) ?
01063                     (rightmost->p.x - tailp->p.x) / (rightmost->p.y - tailp->p.x) : 0; //second right step
01064 
01065     double lstep = flstep, rstep = frstep;
01066 
01067     double left_x = min_y->p.x, right_x = min_y->p.x;
01068 
01069     // Loop around all points in the region and count those that are aligned.
01070     int min_iter = min_y->p.y;
01071     int max_iter = max_y->p.y;
01072     for(int y = min_iter; y <= max_iter; ++y)
01073     {
01074         if (y < 0 || y >= img_height) continue;
01075 
01076         int adx = y * img_width + int(left_x);
01077         for(int x = int(left_x); x <= int(right_x); ++x, ++adx)
01078         {
01079             if (x < 0 || x >= img_width) continue;
01080 
01081             ++total_pts;
01082             if(isAligned(adx, rec.theta, rec.prec))
01083             {
01084                 ++alg_pts;
01085             }
01086         }
01087 
01088         if(y >= leftmost->p.y) { lstep = slstep; }
01089         if(y >= rightmost->p.y) { rstep = srstep; }
01090 
01091         left_x += lstep;
01092         right_x += rstep;
01093     }
01094 
01095     return nfa(total_pts, alg_pts, rec.p);
01096 }
01097 
01098 double LineSegmentDetectorImpl::nfa(const int& n, const int& k, const double& p) const
01099 {
01100     // Trivial cases
01101     if(n == 0 || k == 0) { return -LOG_NT; }
01102     if(n == k) { return -LOG_NT - double(n) * log10(p); }
01103 
01104     double p_term = p / (1 - p);
01105 
01106     double log1term = (double(n) + 1) - log_gamma(double(k) + 1)
01107                 - log_gamma(double(n-k) + 1)
01108                 + double(k) * log(p) + double(n-k) * log(1.0 - p);
01109     double term = exp(log1term);
01110 
01111     if(double_equal(term, 0))
01112     {
01113         if(k > n * p) return -log1term / M_LN10 - LOG_NT;
01114         else return -LOG_NT;
01115     }
01116 
01117     // Compute more terms if needed
01118     double bin_tail = term;
01119     double tolerance = 0.1; // an error of 10% in the result is accepted
01120     for(int i = k + 1; i <= n; ++i)
01121     {
01122         double bin_term = double(n - i + 1) / double(i);
01123         double mult_term = bin_term * p_term;
01124         term *= mult_term;
01125         bin_tail += term;
01126         if(bin_term < 1)
01127         {
01128             double err = term * ((1 - pow(mult_term, double(n-i+1))) / (1 - mult_term) - 1);
01129             if(err < tolerance * fabs(-log10(bin_tail) - LOG_NT) * bin_tail) break;
01130         }
01131 
01132     }
01133     return -log10(bin_tail) - LOG_NT;
01134 }
01135 
01136 inline bool LineSegmentDetectorImpl::isAligned(const int& address, const double& theta, const double& prec) const
01137 {
01138     if(address < 0) { return false; }
01139     const double& a = angles_data[address];
01140     if(a == NOTDEF) { return false; }
01141 
01142     // It is assumed that 'theta' and 'a' are in the range [-pi,pi]
01143     double n_theta = theta - a;
01144     if(n_theta < 0) { n_theta = -n_theta; }
01145     if(n_theta > M_3_2_PI)
01146     {
01147         n_theta -= M_2__PI;
01148         if(n_theta < 0) n_theta = -n_theta;
01149     }
01150 
01151     return n_theta <= prec;
01152 }
01153 
01154 
01155 void LineSegmentDetectorImpl::drawSegments(InputOutputArray _image, InputArray lines)
01156 {
01157     CV_Assert(!_image.empty() && (_image.channels() == 1 || _image.channels() == 3));
01158 
01159     Mat gray;
01160     if (_image.channels() == 1)
01161     {
01162         gray = _image.getMatRef();
01163     }
01164     else if (_image.channels() == 3)
01165     {
01166         cvtColor(_image, gray, CV_BGR2GRAY);
01167     }
01168 
01169     // Create a 3 channel image in order to draw colored lines
01170     std::vector<Mat> planes;
01171     planes.push_back(gray);
01172     planes.push_back(gray);
01173     planes.push_back(gray);
01174 
01175     merge(planes, _image);
01176 
01177     Mat _lines;
01178     _lines = lines.getMat();
01179     int N = _lines.checkVector(4);
01180 
01181     // Draw segments
01182     for(int i = 0; i < N; ++i)
01183     {
01184         const Vec4f& v = _lines.at<Vec4f>(i);
01185         Point2f b(v[0], v[1]);
01186         Point2f e(v[2], v[3]);
01187         line(_image.getMatRef(), b, e, Scalar(0, 0, 255), 1);
01188     }
01189 }
01190 
01191 
01192 int LineSegmentDetectorImpl::compareSegments(const Size& size, InputArray lines1, InputArray lines2, InputOutputArray _image)
01193 {
01194     Size sz = size;
01195     if (_image.needed() && _image.size() != size) sz = _image.size();
01196     CV_Assert(sz.area());
01197 
01198     Mat_<uchar> I1 = Mat_<uchar>::zeros(sz);
01199     Mat_<uchar> I2 = Mat_<uchar>::zeros(sz);
01200 
01201     Mat _lines1;
01202     Mat _lines2;
01203     _lines1 = lines1.getMat();
01204     _lines2 = lines2.getMat();
01205     int N1 = _lines1.checkVector(4);
01206     int N2 = _lines2.checkVector(4);
01207 
01208     // Draw segments
01209     for(int i = 0; i < N1; ++i)
01210     {
01211         Point2f b(_lines1.at<Vec4f>(i)[0], _lines1.at<Vec4f>(i)[1]);
01212         Point2f e(_lines1.at<Vec4f>(i)[2], _lines1.at<Vec4f>(i)[3]);
01213         line(I1, b, e, Scalar::all(255), 1);
01214     }
01215     for(int i = 0; i < N2; ++i)
01216     {
01217         Point2f b(_lines2.at<Vec4f>(i)[0], _lines2.at<Vec4f>(i)[1]);
01218         Point2f e(_lines2.at<Vec4f>(i)[2], _lines2.at<Vec4f>(i)[3]);
01219         line(I2, b, e, Scalar::all(255), 1);
01220     }
01221 
01222     // Count the pixels that don't agree
01223     Mat Ixor;
01224     bitwise_xor(I1, I2, Ixor);
01225     int N = countNonZero(Ixor);
01226 
01227     if (_image.needed())
01228     {
01229         CV_Assert(_image.channels() == 3);
01230         Mat img = _image.getMatRef();
01231         CV_Assert(img.isContinuous() && I1.isContinuous() && I2.isContinuous());
01232 
01233         for (unsigned int i = 0; i < I1.total(); ++i)
01234         {
01235             uchar i1 = I1.ptr()[i];
01236             uchar i2 = I2.ptr()[i];
01237             if (i1 || i2)
01238             {
01239                 unsigned int base_idx = i * 3;
01240                 if (i1) img.ptr()[base_idx] = 255;
01241                 else img.ptr()[base_idx] = 0;
01242                 img.ptr()[base_idx + 1] = 0;
01243                 if (i2) img.ptr()[base_idx + 2] = 255;
01244                 else img.ptr()[base_idx + 2] = 0;
01245             }
01246         }
01247     }
01248 
01249     return N;
01250 }
01251 
01252 } // namespace cv
01253