Merge pull request #1147 from lhelontra:edgenms

Author: vpisarev

modules/ximgproc/include/opencv2/ximgproc/structured_edge_detection.hpp

@@ -102,11 +102,31 @@ class CV_EXPORTS_W StructuredEdgeDetection : public Algorithm
 
     The algorithm underlies this function is much more robust to texture presence, than common
     approaches, e.g. Sobel
-    @param src source image (RGB, float, in [0;1]) to detect edges
-    @param dst destination image (grayscale, float, in [0;1]) where edges are drawn
+    @param _src source image (RGB, float, in [0;1]) to detect edges
+    @param _dst destination image (grayscale, float, in [0;1]) where edges are drawn
     @sa Sobel, Canny
      */
-    CV_WRAP virtual void detectEdges(const Mat &src, CV_OUT Mat &dst) const = 0;
+    CV_WRAP virtual void detectEdges(cv::InputArray _src, cv::OutputArray _dst) const = 0;
+
+    /** @brief The function computes orientation from edge image.
+
+    @param _src edge image.
+    @param _dst orientation image.
+     */
+    CV_WRAP virtual void computeOrientation(cv::InputArray _src, cv::OutputArray _dst) const = 0;
+
+
+    /** @brief The function edgenms in edge image and suppress edges where edge is stronger in orthogonal direction.
+
+    @param edge_image edge image from detectEdges function.
+    @param orientation_image orientation image from computeOrientation function.
+    @param _dst suppressed image (grayscale, float, in [0;1])
+    @param r radius for NMS suppression.
+    @param s radius for boundary suppression.
+    @param m multiplier for conservative suppression.
+    @param isParallel enables/disables parallel computing.
+     */
+    CV_WRAP virtual void edgesNms(cv::InputArray edge_image, cv::InputArray orientation_image, cv::OutputArray _dst, int r = 2, int s = 0, float m = 1, bool isParallel = true) const = 0;
 };
 
 /*!

modules/ximgproc/samples/structured_edge_detection.cpp

@@ -59,10 +59,20 @@ int main( int argc, const char** argv )
         createStructuredEdgeDetection(modelFilename);
     pDollar->detectEdges(image, edges);
 
+    // computes orientation from edge map
+    Mat orientation_map;
+    pDollar->computeOrientation(edges, orientation_map);
+
+    // suppress edges
+    Mat edge_nms;
+    pDollar->edgesNms(edges, orientation_map, edge_nms, 2, 0, 1, true);
+
     if ( outFilename.size() == 0 )
     {
         namedWindow("edges", 1);
         imshow("edges", edges);
+        namedWindow("edges nms", 1);
+        imshow("edges nms", edge_nms);
         waitKey(0);
     }
     else

modules/ximgproc/src/structured_edge_detection.cpp

@@ -255,6 +255,81 @@ static void gradientHist(const cv::Mat &src, cv::Mat &magnitude, cv::Mat &histog
     }
 }
 
+/*!
+ * The class parallelizing the edgenms algorithm.
+ *
+ * \param E : edge image
+ * \param O : orientation image
+ * \param dst : destination image
+ * \param r : radius for NMS suppression
+ * \param s : radius for boundary suppression
+ * \param m : multiplier for conservative suppression
+ */
+class NmsInvoker : public cv::ParallelLoopBody
+{
+
+private:
+  const cv::Mat &E;
+  const cv::Mat &O;
+  cv::Mat &dst;
+  const int r;
+  const float m;
+
+public:
+  NmsInvoker(const cv::Mat &_E, const cv::Mat &_O, cv::Mat &_dst, const int _r, const float _m)
+              : E(_E), O(_O), dst(_dst), r(_r), m(_m)
+              {
+              }
+
+  void operator()(const cv::Range &range) const
+  {
+     for (int x = range.start; x < range.end; x++)
+     {
+       const float *e_ptr = E.ptr<float>(x);
+       const float *o_ptr = O.ptr<float>(x);
+       float *dst_ptr = dst.ptr<float>(x);
+       for (int y=0; y < E.cols; y++)
+       {
+         float e = e_ptr[y];
+         dst_ptr[y] = e;
+         if (!e) continue;
+         e *= m;
+         float coso = cos(o_ptr[y]);
+         float sino = sin(o_ptr[y]);
+         for (int d=-r; d<=r; d++)
+         {
+           if (d)
+           {
+             float xdcos = x+d*coso;
+             float ydsin = y+d*sino;
+             xdcos = xdcos < 0 ? 0 : (xdcos > E.rows - 1.001f ? E.rows - 1.001f : xdcos);
+             ydsin = ydsin < 0 ? 0 : (ydsin > E.cols - 1.001f ? E.cols - 1.001f : ydsin);
+             int x0 = (int)xdcos;
+             int y0 = (int)ydsin;
+             int x1 = x0 + 1;
+             int y1 = y0 + 1;
+             float dx0 = xdcos - x0;
+             float dy0 = ydsin - y0;
+             float dx1 = 1 - dx0;
+             float dy1 = 1 - dy0;
+             float e0 = E.at<float>(x0, y0) * dx1 * dy1 +
+                         E.at<float>(x1, y0) * dx0 * dy1 +
+                         E.at<float>(x0, y1) * dx1 * dy0 +
+                         E.at<float>(x1, y1) * dx0 * dy0;
+
+             if(e < e0)
+             {
+               dst_ptr[y] = 0;
+               break;
+             }
+           }
+         }
+
+       }
+     }
+  }
+};
+
 /********************* RFFeatureGetter class *********************/
 
 namespace cv
@@ -445,21 +520,23 @@ class StructuredEdgeDetectionImpl : public StructuredEdgeDetection
     /*!
      * The function detects edges in src and draw them to dst
      *
-     * \param src : source image (RGB, float, in [0;1]) to detect edges
-     * \param dst : destination image (grayscale, float, in [0;1])
+     * \param _src : source image (RGB, float, in [0;1]) to detect edges
+     * \param _dst : destination image (grayscale, float, in [0;1])
      *              where edges are drawn
      */
-    void detectEdges(const cv::Mat &src, cv::Mat &dst) const
+    void detectEdges(cv::InputArray _src, cv::OutputArray _dst) const
     {
-        CV_Assert( src.type() == CV_32FC3 );
+        CV_Assert( _src.type() == CV_32FC3 );
 
-        dst.create( src.size(), cv::DataType<float>::type );
+        _dst.createSameSize( _src, cv::DataType<float>::type );
+        _dst.setTo(0);
+        Mat dst = _dst.getMat();
 
         int padding = ( __rf.options.patchSize
             - __rf.options.patchInnerSize )/2;
 
         cv::Mat nSrc;
-        copyMakeBorder( src, nSrc, padding, padding,
+        copyMakeBorder( _src, nSrc, padding, padding,
             padding, padding, BORDER_REFLECT );
 
         NChannelsMat features;
@@ -472,6 +549,101 @@ class StructuredEdgeDetectionImpl : public StructuredEdgeDetection
         predictEdges( features, dst );
     }
 
+    /*!
+     * The function computes orientation from edge image.
+     *
+     * \param src : edge image.
+     * \param dst : orientation image.
+     * \param r : filter radius.
+     */
+    void computeOrientation(cv::InputArray _src, cv::OutputArray _dst) const
+    {
+      CV_Assert( _src.type() == CV_32FC1 );
+
+      cv::Mat Oxx, Oxy, Oyy;
+
+      _dst.createSameSize( _src, _src.type() );
+      _dst.setTo(0);
+
+      Mat src = _src.getMat();
+      cv::Mat E_conv = imsmooth(src, __rf.options.gradientNormalizationRadius);
+
+      Sobel(E_conv, Oxx, -1, 2, 0);
+      Sobel(E_conv, Oxy, -1, 1, 1);
+      Sobel(E_conv, Oyy, -1, 0, 2);
+
+      Mat dst = _dst.getMat();
+      float *o = dst.ptr<float>();
+      float *oxx = Oxx.ptr<float>();
+      float *oxy = Oxy.ptr<float>();
+      float *oyy = Oyy.ptr<float>();
+      for (int i = 0; i < dst.rows * dst.cols; i++)
+      {
+          int xysign = -((oxy[i] > 0) - (oxy[i] < 0));
+          o[i] = (atan((oyy[i] * xysign / (oxx[i] + 1e-5))) > 0) ? (float) fmod(
+                  atan((oyy[i] * xysign / (oxx[i] + 1e-5))), M_PI) : (float) fmod(
+                  atan((oyy[i] * xysign / (oxx[i] + 1e-5))) + M_PI, M_PI);
+      }
+    }
+
+     /*!
+     * The function suppress edges where edge is stronger in orthogonal direction
+     * \param edge_image : edge image from detectEdges function.
+     * \param orientation_image : orientation image from computeOrientation function.
+     * \param _dst : suppressed image (grayscale, float, in [0;1])
+     * \param r : radius for NMS suppression.
+     * \param s : radius for boundary suppression.
+     * \param m : multiplier for conservative suppression.
+     * \param isParallel: enables/disables parallel computing.
+     */
+    void edgesNms(cv::InputArray edge_image, cv::InputArray orientation_image, cv::OutputArray _dst, int r, int s, float m, bool isParallel) const
+    {
+        CV_Assert(edge_image.type() == CV_32FC1);
+        CV_Assert(orientation_image.type() == CV_32FC1);
+
+        cv::Mat E = edge_image.getMat();
+        cv::Mat O = orientation_image.getMat();
+        cv::Mat E_t = E.t();
+        cv::Mat O_t = O.t();
+
+        cv::Mat dst = _dst.getMat();
+        dst.create(E.cols, E.rows, E.type());
+        dst.setTo(0);
+
+        cv::Range sizeRange = cv::Range(0, E_t.rows);
+        NmsInvoker body = NmsInvoker(E_t, O_t, dst, r, m);
+        if (isParallel)
+        {
+          cv::parallel_for_(sizeRange, body);
+        } else
+        {
+          body(sizeRange);
+        }
+
+        s = s > E_t.rows / 2 ? E_t.rows / 2 : s;
+        s = s > E_t.cols / 2 ? E_t.cols / 2 : s;
+        for (int x=0; x<s; x++)
+        {
+          for (int y=0; y<E_t.cols; y++)
+          {
+            dst.at<float>(x, y) *= x / (float)s;
+            dst.at<float>(E_t.rows-1-x, y) *= x / (float)s;
+          }
+        }
+
+        for (int x=0; x < E_t.rows; x++)
+        {
+          for (int y=0; y < s; y++)
+          {
+            dst.at<float>(x, y) *= y / (float)s;
+            dst.at<float>(x, E_t.cols-1-y) *= y / (float)s;
+          }
+        }
+      transpose(dst, dst);
+      dst.copyTo(_dst);
+    }
+
+
 protected:
     /*!
      * Private method used by process method. The function