Merge pull request #7311 from fran6co:lsd

modules/imgproc/src/lsd.cpp

@@ -230,12 +230,9 @@ public:
 
 private:
     Mat image;
-    Mat_<double> scaled_image;
-    double *scaled_image_data;
+    Mat scaled_image;
     Mat_<double> angles;     // in rads
-    double *angles_data;
     Mat_<double> modgrad;
-    double *modgrad_data;
     Mat_<uchar> used;
 
     int img_width;
@@ -270,6 +267,8 @@ private:
         struct coorlist* next;
     };
 
+    std::vector<coorlist> list;
+
     struct rect
     {
         double x1, y1, x2, y2;    // first and second point of the line segment
@@ -309,7 +308,7 @@ private:
  * @param list      Return: Vector of coordinate points that are pseudo ordered by magnitude.
  *                  Pixels would be ordered by norm value, up to a precision given by max_grad/n_bins.
  */
-    void ll_angle(const double& threshold, const unsigned int& n_bins, std::vector<coorlist>& list);
+    void ll_angle(const double& threshold, const unsigned int& n_bins);
 
 /**
  * Grow a region starting from point s with a defined precision,
@@ -317,31 +316,29 @@ private:
  *
  * @param s         Starting point for the region.
  * @param reg       Return: Vector of points, that are part of the region
- * @param reg_size  Return: The size of the region.
  * @param reg_angle Return: The mean angle of the region.
  * @param prec      The precision by which each region angle should be aligned to the mean.
  */
     void region_grow(const Point2i& s, std::vector<RegionPoint>& reg,
-                     int& reg_size, double& reg_angle, const double& prec);
+                     double& reg_angle, const double& prec);
 
 /**
  * Finds the bounding rotated rectangle of a region.
  *
  * @param reg       The region of points, from which the rectangle to be constructed from.
- * @param reg_size  The number of points in the region.
  * @param reg_angle The mean angle of the region.
  * @param prec      The precision by which points were found.
  * @param p         Probability of a point with angle within 'prec'.
  * @param rec       Return: The generated rectangle.
  */
-    void region2rect(const std::vector<RegionPoint>& reg, const int reg_size, const double reg_angle,
+    void region2rect(const std::vector<RegionPoint>& reg, const double reg_angle,
                      const double prec, const double p, rect& rec) const;
 
 /**
  * Compute region's angle as the principal inertia axis of the region.
  * @return          Regions angle.
  */
-    double get_theta(const std::vector<RegionPoint>& reg, const int& reg_size, const double& x,
+    double get_theta(const std::vector<RegionPoint>& reg, const double& x,
                      const double& y, const double& reg_angle, const double& prec) const;
 
 /**
@@ -350,14 +347,14 @@ private:
  * estimated angle tolerance. If this fails to produce a rectangle with the right density of region points,
  * 'reduce_region_radius' is called to try to satisfy this condition.
  */
-    bool refine(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
+    bool refine(std::vector<RegionPoint>& reg, double reg_angle,
                 const double prec, double p, rect& rec, const double& density_th);
 
 /**
  * Reduce the region size, by elimination the points far from the starting point, until that leads to
  * rectangle with the right density of region points or to discard the region if too small.
  */
-    bool reduce_region_radius(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
+    bool reduce_region_radius(std::vector<RegionPoint>& reg, double reg_angle,
                 const double prec, double p, rect& rec, double density, const double& density_th);
 
 /**
@@ -383,7 +380,7 @@ private:
  * Is the point at place 'address' aligned to angle theta, up to precision 'prec'?
  * @return      Whether the point is aligned.
  */
-    bool isAligned(const int& address, const double& theta, const double& prec) const;
+    bool isAligned(int x, int y, const double& theta, const double& prec) const;
 };
 
 /////////////////////////////////////////////////////////////////////////////////////////
@@ -414,11 +411,8 @@ void LineSegmentDetectorImpl::detect(InputArray _image, OutputArray _lines,
 {
     CV_INSTRUMENT_REGION()
 
-    Mat_<double> img = _image.getMat();
-    CV_Assert(!img.empty() && img.channels() == 1);
-
-    // Convert image to double
-    img.convertTo(image, CV_64FC1);
+    image = _image.getMat();
+    CV_Assert(!image.empty() && image.type() == CV_8UC1);
 
     std::vector<Vec4f> lines;
     std::vector<double> w, p, n;
@@ -446,7 +440,6 @@ void LineSegmentDetectorImpl::flsd(std::vector<Vec4f>& lines,
     const double p = ANG_TH / 180;
     const double rho = QUANT / sin(prec);    // gradient magnitude threshold
 
-    std::vector<coorlist> list;
     if(SCALE != 1)
     {
         Mat gaussian_img;
@@ -457,45 +450,43 @@ void LineSegmentDetectorImpl::flsd(std::vector<Vec4f>& lines,
         GaussianBlur(image, gaussian_img, ksize, sigma);
         // Scale image to needed size
         resize(gaussian_img, scaled_image, Size(), SCALE, SCALE);
-        ll_angle(rho, N_BINS, list);
+        ll_angle(rho, N_BINS);
     }
     else
     {
         scaled_image = image;
-        ll_angle(rho, N_BINS, list);
+        ll_angle(rho, N_BINS);
     }
 
     LOG_NT = 5 * (log10(double(img_width)) + log10(double(img_height))) / 2 + log10(11.0);
-    const int min_reg_size = int(-LOG_NT/log10(p)); // minimal number of points in region that can give a meaningful event
+    const size_t min_reg_size = size_t(-LOG_NT/log10(p)); // minimal number of points in region that can give a meaningful event
 
     // // Initialize region only when needed
     // Mat region = Mat::zeros(scaled_image.size(), CV_8UC1);
     used = Mat_<uchar>::zeros(scaled_image.size()); // zeros = NOTUSED
-    std::vector<RegionPoint> reg(img_width * img_height);
+    std::vector<RegionPoint> reg;
 
     // Search for line segments
-    unsigned int ls_count = 0;
     for(size_t i = 0, list_size = list.size(); i < list_size; ++i)
     {
-        unsigned int adx = list[i].p.x + list[i].p.y * img_width;
-        if((used.ptr()[adx] == NOTUSED) && (angles_data[adx] != NOTDEF))
+        const Point2i& point = list[i].p;
+        if((used.at<uchar>(point) == NOTUSED) && (angles.at<double>(point) != NOTDEF))
         {
-            int reg_size;
             double reg_angle;
-            region_grow(list[i].p, reg, reg_size, reg_angle, prec);
+            region_grow(list[i].p, reg, reg_angle, prec);
 
             // Ignore small regions
-            if(reg_size < min_reg_size) { continue; }
+            if(reg.size() < min_reg_size) { continue; }
 
             // Construct rectangular approximation for the region
             rect rec;
-            region2rect(reg, reg_size, reg_angle, prec, p, rec);
+            region2rect(reg, reg_angle, prec, p, rec);
 
             double log_nfa = -1;
             if(doRefine > LSD_REFINE_NONE)
             {
                 // At least REFINE_STANDARD lvl.
-                if(!refine(reg, reg_size, reg_angle, prec, p, rec, DENSITY_TH)) { continue; }
+                if(!refine(reg, reg_angle, prec, p, rec, DENSITY_TH)) { continue; }
 
                 if(doRefine >= LSD_REFINE_ADV)
                 {
@@ -505,7 +496,6 @@ void LineSegmentDetectorImpl::flsd(std::vector<Vec4f>& lines,
                 }
             }
             // Found new line
-            ++ls_count;
 
             // Add the offset
             rec.x1 += 0.5; rec.y1 += 0.5;
@@ -524,29 +514,17 @@ void LineSegmentDetectorImpl::flsd(std::vector<Vec4f>& lines,
             if(w_needed) widths.push_back(rec.width);
             if(p_needed) precisions.push_back(rec.p);
             if(n_needed && doRefine >= LSD_REFINE_ADV) nfas.push_back(log_nfa);
-
-
-            // //Add the linesID to the region on the image
-            // for(unsigned int el = 0; el < reg_size; el++)
-            // {
-            //     region.data[reg[i].x + reg[i].y * width] = ls_count;
-            // }
         }
     }
 }
 
 void LineSegmentDetectorImpl::ll_angle(const double& threshold,
-                                   const unsigned int& n_bins,
-                                   std::vector<coorlist>& list)
+                                   const unsigned int& n_bins)
 {
     //Initialize data
     angles = Mat_<double>(scaled_image.size());
     modgrad = Mat_<double>(scaled_image.size());
 
-    angles_data = angles.ptr<double>(0);
-    modgrad_data = modgrad.ptr<double>(0);
-    scaled_image_data = scaled_image.ptr<double>(0);
-
     img_width = scaled_image.cols;
     img_height = scaled_image.rows;
 
@@ -555,30 +533,30 @@ void LineSegmentDetectorImpl::ll_angle(const double& threshold,
     angles.col(img_width - 1).setTo(NOTDEF);
 
     // Computing gradient for remaining pixels
-    CV_Assert(scaled_image.isContinuous() &&
-              modgrad.isContinuous() &&
-              angles.isContinuous());   // Accessing image data linearly
-
     double max_grad = -1;
     for(int y = 0; y < img_height - 1; ++y)
     {
-        for(int addr = y * img_width, addr_end = addr + img_width - 1; addr < addr_end; ++addr)
+        const uchar* scaled_image_row = scaled_image.ptr<uchar>(y);
+        const uchar* next_scaled_image_row = scaled_image.ptr<uchar>(y+1);
+        double* angles_row = angles.ptr<double>(y);
+        double* modgrad_row = modgrad.ptr<double>(y);
+        for(int x = 0; x < img_width-1; ++x)
         {
-            double DA = scaled_image_data[addr + img_width + 1] - scaled_image_data[addr];
-            double BC = scaled_image_data[addr + 1] - scaled_image_data[addr + img_width];
-            double gx = DA + BC;    // gradient x component
-            double gy = DA - BC;    // gradient y component
-            double norm = std::sqrt((gx * gx + gy * gy) / 4); // gradient norm
+            int DA = next_scaled_image_row[x + 1] - scaled_image_row[x];
+            int BC = scaled_image_row[x + 1] - next_scaled_image_row[x];
+            int gx = DA + BC;    // gradient x component
+            int gy = DA - BC;    // gradient y component
+            double norm = std::sqrt((gx * gx + gy * gy) / 4.0); // gradient norm
 
-            modgrad_data[addr] = norm;    // store gradient
+            modgrad_row[x] = norm;    // store gradient
 
             if (norm <= threshold)  // norm too small, gradient no defined
             {
-                angles_data[addr] = NOTDEF;
+                angles_row[x] = NOTDEF;
             }
             else
             {
-                angles_data[addr] = fastAtan2(float(gx), float(-gy)) * DEG_TO_RADS;  // gradient angle computation
+                angles_row[x] = fastAtan2(float(gx), float(-gy)) * DEG_TO_RADS;  // gradient angle computation
                 if (norm > max_grad) { max_grad = norm; }
             }
 
@@ -586,7 +564,7 @@ void LineSegmentDetectorImpl::ll_angle(const double& threshold,
     }
 
     // Compute histogram of gradient values
-    list = std::vector<coorlist>(img_width * img_height);
+    list.resize(img_width * img_height);
     std::vector<coorlist*> range_s(n_bins);
     std::vector<coorlist*> range_e(n_bins);
     unsigned int count = 0;
@@ -594,11 +572,11 @@ void LineSegmentDetectorImpl::ll_angle(const double& threshold,
 
     for(int y = 0; y < img_height - 1; ++y)
     {
-        const double* norm = modgrad_data + y * img_width;
-        for(int x = 0; x < img_width - 1; ++x, ++norm)
+        const double* modgrad_row = modgrad.ptr<double>(y);
+        for(int x = 0; x < img_width - 1; ++x)
         {
             // Store the point in the right bin according to its norm
-            int i = int((*norm) * bin_coef);
+            int i = int(modgrad_row[x] * bin_coef);
             if(!range_e[i])
             {
                 range_e[i] = range_s[i] = &list[count];
@@ -635,46 +613,51 @@ void LineSegmentDetectorImpl::ll_angle(const double& threshold,
 }
 
 void LineSegmentDetectorImpl::region_grow(const Point2i& s, std::vector<RegionPoint>& reg,
-                                      int& reg_size, double& reg_angle, const double& prec)
+                                      double& reg_angle, const double& prec)
 {
+    reg.clear();
+
     // Point to this region
-    reg_size = 1;
-    reg[0].x = s.x;
-    reg[0].y = s.y;
-    int addr = s.x + s.y * img_width;
-    reg[0].used = used.ptr() + addr;
-    reg_angle = angles_data[addr];
-    reg[0].angle = reg_angle;
-    reg[0].modgrad = modgrad_data[addr];
+    RegionPoint seed;
+    seed.x = s.x;
+    seed.y = s.y;
+    seed.used = &used.at<uchar>(s);
+    reg_angle = angles.at<double>(s);
+    seed.angle = reg_angle;
+    seed.modgrad = modgrad.at<double>(s);
+    reg.push_back(seed);
 
     float sumdx = float(std::cos(reg_angle));
     float sumdy = float(std::sin(reg_angle));
-    *reg[0].used = USED;
+    *seed.used = USED;
 
     //Try neighboring regions
-    for(int i = 0; i < reg_size; ++i)
+    for (size_t i = 0;i<reg.size();i++)
     {
         const RegionPoint& rpoint = reg[i];
         int xx_min = std::max(rpoint.x - 1, 0), xx_max = std::min(rpoint.x + 1, img_width - 1);
         int yy_min = std::max(rpoint.y - 1, 0), yy_max = std::min(rpoint.y + 1, img_height - 1);
         for(int yy = yy_min; yy <= yy_max; ++yy)
         {
-            int c_addr = xx_min + yy * img_width;
-            for(int xx = xx_min; xx <= xx_max; ++xx, ++c_addr)
+            uchar* used_row = used.ptr<uchar>(yy);
+            const double* angles_row = angles.ptr<double>(yy);
+            const double* modgrad_row = modgrad.ptr<double>(yy);
+            for(int xx = xx_min; xx <= xx_max; ++xx)
             {
-                if((used.ptr()[c_addr] != USED) &&
-                   (isAligned(c_addr, reg_angle, prec)))
+                uchar& is_used = used_row[xx];
+                if(is_used != USED &&
+                   (isAligned(xx, yy, reg_angle, prec)))
                 {
+                    const double& angle = angles_row[xx];
                     // Add point
-                    used.ptr()[c_addr] = USED;
-                    RegionPoint& region_point = reg[reg_size];
+                    is_used = USED;
+                    RegionPoint region_point;
                     region_point.x = xx;
                     region_point.y = yy;
-                    region_point.used = &(used.ptr()[c_addr]);
-                    region_point.modgrad = modgrad_data[c_addr];
-                    const double& angle = angles_data[c_addr];
+                    region_point.used = &is_used;
+                    region_point.modgrad = modgrad_row[xx];
                     region_point.angle = angle;
-                    ++reg_size;
+                    reg.push_back(region_point);
 
                     // Update region's angle
                     sumdx += cos(float(angle));
@@ -687,11 +670,11 @@ void LineSegmentDetectorImpl::region_grow(const Point2i& s, std::vector<RegionPo
     }
 }
 
-void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg, const int reg_size,
+void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg,
                                       const double reg_angle, const double prec, const double p, rect& rec) const
 {
     double x = 0, y = 0, sum = 0;
-    for(int i = 0; i < reg_size; ++i)
+    for(size_t i = 0; i < reg.size(); ++i)
     {
         const RegionPoint& pnt = reg[i];
         const double& weight = pnt.modgrad;
@@ -706,14 +689,14 @@ void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg, c
     x /= sum;
     y /= sum;
 
-    double theta = get_theta(reg, reg_size, x, y, reg_angle, prec);
+    double theta = get_theta(reg, x, y, reg_angle, prec);
 
     // Find length and width
     double dx = cos(theta);
     double dy = sin(theta);
     double l_min = 0, l_max = 0, w_min = 0, w_max = 0;
 
-    for(int i = 0; i < reg_size; ++i)
+    for(size_t i = 0; i < reg.size(); ++i)
     {
         double regdx = double(reg[i].x) - x;
         double regdy = double(reg[i].y) - y;
@@ -745,7 +728,7 @@ void LineSegmentDetectorImpl::region2rect(const std::vector<RegionPoint>& reg, c
     if(rec.width < 1.0) rec.width = 1.0;
 }
 
-double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, const int& reg_size, const double& x,
+double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, const double& x,
                                       const double& y, const double& reg_angle, const double& prec) const
 {
     double Ixx = 0.0;
@@ -753,7 +736,7 @@ double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, c
     double Ixy = 0.0;
 
     // Compute inertia matrix
-    for(int i = 0; i < reg_size; ++i)
+    for(size_t i = 0; i < reg.size(); ++i)
     {
         const double& regx = reg[i].x;
         const double& regy = reg[i].y;
@@ -783,10 +766,10 @@ double LineSegmentDetectorImpl::get_theta(const std::vector<RegionPoint>& reg, c
     return theta;
 }
 
-bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
+bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, double reg_angle,
                                  const double prec, double p, rect& rec, const double& density_th)
 {
-    double density = double(reg_size) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
+    double density = double(reg.size()) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
 
     if (density >= density_th) { return true; }
 
@@ -797,7 +780,7 @@ bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, int& reg_siz
     double sum = 0, s_sum = 0;
     int n = 0;
 
-    for (int i = 0; i < reg_size; ++i)
+    for (size_t i = 0; i < reg.size(); ++i)
     {
         *(reg[i].used) = NOTUSED;
         if (dist(xc, yc, reg[i].x, reg[i].y) < rec.width)
@@ -814,16 +797,16 @@ bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, int& reg_siz
     double tau = 2.0 * sqrt((s_sum - 2.0 * mean_angle * sum) / double(n) + mean_angle * mean_angle);
 
     // Try new region
-    region_grow(Point(reg[0].x, reg[0].y), reg, reg_size, reg_angle, tau);
+    region_grow(Point(reg[0].x, reg[0].y), reg, reg_angle, tau);
 
-    if (reg_size < 2) { return false; }
+    if (reg.size() < 2) { return false; }
 
-    region2rect(reg, reg_size, reg_angle, prec, p, rec);
-    density = double(reg_size) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
+    region2rect(reg, reg_angle, prec, p, rec);
+    density = double(reg.size()) / (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
 
     if (density < density_th)
     {
-        return reduce_region_radius(reg, reg_size, reg_angle, prec, p, rec, density, density_th);
+        return reduce_region_radius(reg, reg_angle, prec, p, rec, density, density_th);
     }
     else
     {
@@ -831,7 +814,7 @@ bool LineSegmentDetectorImpl::refine(std::vector<RegionPoint>& reg, int& reg_siz
     }
 }
 
-bool LineSegmentDetectorImpl::reduce_region_radius(std::vector<RegionPoint>& reg, int& reg_size, double reg_angle,
+bool LineSegmentDetectorImpl::reduce_region_radius(std::vector<RegionPoint>& reg, double reg_angle,
                 const double prec, double p, rect& rec, double density, const double& density_th)
 {
     // Compute region's radius
@@ -845,25 +828,25 @@ bool LineSegmentDetectorImpl::reduce_region_radius(std::vector<RegionPoint>& reg
     {
         radSq *= 0.75*0.75; // Reduce region's radius to 75% of its value
         // Remove points from the region and update 'used' map
-        for(int i = 0; i < reg_size; ++i)
+        for (size_t i = 0; i < reg.size(); ++i)
         {
             if(distSq(xc, yc, double(reg[i].x), double(reg[i].y)) > radSq)
             {
                 // Remove point from the region
                 *(reg[i].used) = NOTUSED;
-                std::swap(reg[i], reg[reg_size - 1]);
-                --reg_size;
+                std::swap(reg[i], reg[reg.size() - 1]);
+                reg.pop_back();
                 --i; // To avoid skipping one point
             }
         }
 
-        if(reg_size < 2) { return false; }
+        if(reg.size() < 2) { return false; }
 
         // Re-compute rectangle
-        region2rect(reg, reg_size ,reg_angle, prec, p, rec);
+        region2rect(reg ,reg_angle, prec, p, rec);
 
         // Re-compute region points density
-        density = double(reg_size) /
+        density = double(reg.size()) /
                   (dist(rec.x1, rec.y1, rec.x2, rec.y2) * rec.width);
     }
 
@@ -981,7 +964,7 @@ double LineSegmentDetectorImpl::rect_nfa(const rect& rec) const
     double dyhw = rec.dy * half_width;
     double dxhw = rec.dx * half_width;
 
-    std::vector<edge> ordered_x(4);
+    edge ordered_x[4];
     edge* min_y = &ordered_x[0];
     edge* max_y = &ordered_x[0]; // Will be used for loop range
 
@@ -990,7 +973,7 @@ double LineSegmentDetectorImpl::rect_nfa(const rect& rec) const
     ordered_x[2].p.x = int(rec.x2 + dyhw); ordered_x[2].p.y = int(rec.y2 - dxhw); ordered_x[2].taken = false;
     ordered_x[3].p.x = int(rec.x1 + dyhw); ordered_x[3].p.y = int(rec.y1 - dxhw); ordered_x[3].taken = false;
 
-    std::sort(ordered_x.begin(), ordered_x.end(), AsmallerB_XoverY);
+    std::sort(ordered_x, ordered_x + 4, AsmallerB_XoverY);
 
     // Find min y. And mark as taken. find max y.
     for(unsigned int i = 1; i < 4; ++i)
@@ -1075,13 +1058,12 @@ double LineSegmentDetectorImpl::rect_nfa(const rect& rec) const
     {
         if (y < 0 || y >= img_height) continue;
 
-        int adx = y * img_width + int(left_x);
-        for(int x = int(left_x); x <= int(right_x); ++x, ++adx)
+        for(int x = int(left_x); x <= int(right_x); ++x)
         {
             if (x < 0 || x >= img_width) continue;
 
             ++total_pts;
-            if(isAligned(adx, rec.theta, rec.prec))
+            if(isAligned(x, y, rec.theta, rec.prec))
             {
                 ++alg_pts;
             }
@@ -1135,10 +1117,10 @@ double LineSegmentDetectorImpl::nfa(const int& n, const int& k, const double& p)
     return -log10(bin_tail) - LOG_NT;
 }
 
-inline bool LineSegmentDetectorImpl::isAligned(const int& address, const double& theta, const double& prec) const
+inline bool LineSegmentDetectorImpl::isAligned(int x, int y, const double& theta, const double& prec) const
 {
-    if(address < 0) { return false; }
-    const double& a = angles_data[address];
+    if(x < 0 || y < 0 || x >= angles.cols || y >= angles.rows) { return false; }
+    const double& a = angles.at<double>(y, x);
     if(a == NOTDEF) { return false; }
 
     // It is assumed that 'theta' and 'a' are in the range [-pi,pi]