Skip to content

D
DetectCircle
Commit 3676ea7d authored by 王芳良's avatar 王芳良
Browse files
Options

init

parents
Pipeline #874 canceled with stages
      Show whitespace changes
      Inline Side-by-side
      Showing with 910 additions and 0 deletions
      .gitignore 0 → 100644
      View file @ 3676ea7d
      build/
      bin/
      .vscode/
      .VSCodeCounter/
      *.DS_Store
      *.user
      *.autosave
      *.o
      *.doc#
      Library/
      CMakeLists.txt 0 → 100644
      View file @ 3676ea7d
      cmake_minimum_required(VERSION 3.16)
      PROJECT (Ellipse)
      set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O0 -g")
      find_package(VTK REQUIRED)
      include(${VTK_USE_FILE})
      # INCLUDE_DIRECTORIES(${VTK_USE_FILE})
      find_package(OpenCV REQUIRED)
      message(STATUS "OpenCV library status:")
      message(STATUS " version: ${OpenCV_VERSION}")
      message(STATUS " libraries: ${OpenCV_LIBS}")
      message(STATUS " include path: ${OpenCV_INCLUDE_DIRS}")
      link_directories(
      ${QT_LIBRARY_DIR})
      find_package(Qt5Core REQUIRED)
      # 2. 此处设置为自己的cpp文件，直接添加即可
      SET(SRC_FILES
      main.cpp
      myHough.cpp
      )
      set(HEAD_FILES
      myHough.h
      )
      add_executable(Ellipse ${SRC_FILES} ${HEAD_FILES})
      target_link_libraries(Ellipse ${OpenCV_LIBS} ${VTK_LIBRARIES} Qt5::Core)
      SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin)
      ReadMe.md 0 → 100644
      View file @ 3676ea7d
      ### 使用方法：
      1.在本地bin文件夹下放置含有dcm文件的文件夹foldName
      2.编译
      mkdir build
      cd build
      cmake ..
      make
      3.执行,检测结果保存在Result文件夹内
      cd ../bin
      mkdir Result
      ./Ellipse foldName
      ### 改进措施
      1.先后高斯滤波(去除噪声)、双边滤波(保边去噪)的图像预处理
      2.修改圆检测最大圆半径为40(某些图像无法检测到大圆)
      3.修改圆检测的源码，将HoughCircleEstimateRadiusInvoker函数中的accThreshold传入参数改为原先的一半：
      效果：累加器->潜在圆心：仍使用accThreshold确保数量可控
      潜在圆心->计算半径后的圆心：使用accThreshold/2,防止某些小圆无法被识别
      main.cpp 0 → 100644
      View file @ 3676ea7d
      #include <opencv2/opencv.hpp>
      #include <vector>
      #include <QStandardPaths>
      #include <QDir>
      #include <QDebug>
      #include "vtkSmartPointer.h"
      #include "vtkDICOMImageReader.h"
      #include "vtkImageCast.h"
      #include "vtkImageData.h"
      #include "myHough.h"
      /*
      标准椭圆：x^2/a^2+y^2/b^=1(a>b>0)---------a为长轴长，b为短轴长
      一般椭圆：(可由标准椭圆先旋转angle,再平移center得到)---angle为长轴相对于X轴的偏移角
      Ax^2+Bxy+Cy^2+Dx+Ey+F=0
      故若center为(0,0),则系数D、E=0；若angle为0，则系数B=0。
      将a,b,center,angle代入从一般椭圆转回标准椭圆的过程，即可得到ABCDE->a,b,center,angle的公式
      https://blog.csdn.net/zh471021698/article/details/108812050
      */
      /*
      opencv支持读取dcm方式：
      1.sudo apt-get install libgdcm2-dev
      2.ccmake中将 WITH_GDCM设置为on
      3.重新编译安装
      */
      struct CircleCenter
      {
      cv::Vec3f circle;
      int distance;
      };
      bool mycompare(CircleCenter a, CircleCenter b)
      {
      return a.distance < b.distance; // 降序排列
      }
      void PreProcessing(const cv::Mat &gray, cv::Mat &blurred)
      {
      /*高斯滤波:https://blog.csdn.net/wuqindeyunque/article/details/103694900
      opencv实现的高斯滤波，是对传入的sigmaX，sigmaY分别产生两个一维卷积核，然后先后在行和列上做卷积;
      ==cv2.sepFilter2D(image_gray,-1,cv2.getGaussianKernel(5,1.5),cv2.getGaussianKernel(5,1.5))
      cv::Mat gaussianKernel = cv::getGaussianKernel(5, 1.5);
      //0.1200783842432135; 0.2338807565853503; 0.2920817183428724; 0.2338807565853503; 0.1200783842432135
      cv::getGaussianKernel(5, 3)
      //0.1782032576265784; 0.2105222740037377; 0.2225489367393678; 0.2105222740037377; 0.1782032576265784
      @cv::Size(5, 5)：
      Gaussian kernel size. ksize.width and ksize.height can differ but they both must be positive and odd.
      @1.5:
      Gaussian kernel standard deviation in X direction.if sigmaY is zero, it is set to beequal to sigmaX.
      */
      cv::Mat gaussianBlured;
      cv::GaussianBlur(gray, gaussianBlured, cv::Size(5, 5), 1.5); // 高斯滤波:高斯核进行卷积，平滑边界，去除噪声
      /*双边滤波
      @5：
      d,在滤波时选取的空间距离参数，这里表示以当前像素点为中心点的直径
      @150:
      sigmaColor,是滤波处理时选取的颜色差值范围，该值决定了周围哪些像素点能够参与到滤波中来。
      与当前像素点的像素值差值小于 sigmaColor 的像素点，能够参与到当前的滤波中。该值越大，就说明周围有越多的像素点可以参与到运算中
      @150:
      sigmaSpace,是坐标空间中的sigma值。
      它的值越大，说明有越多的点能够参与到滤波计算中来。当d>0时，无论sigmaSpace的值如何，d都指定邻域大小；否则，d与 sigmaSpace的值成比例。
      */
      cv::Mat bilateralBlured;
      cv::bilateralFilter(gaussianBlured, bilateralBlured, 5, 150, 150);
      blurred = bilateralBlured.clone();
      }
      int HoughCirclesExample(std::string fileName, std::string foldName, std::string fileType)
      {
      // 步骤1：读取图像
      // cv::Mat src = cv::imread(argv[1]);
      cv::Mat src = cv::imread(foldName + "/" + fileName + fileType, cv::IMREAD_GRAYSCALE); // src.type()== 2 == CV_16UC1
      if (src.empty())
      {
      std::cout << "Could not read the image" << std::endl;
      return 1;
      }
      cv::normalize(src, src, 0, 255, cv::NORM_MINMAX, CV_8UC1); // src.type()== 0 == CV_8UC1
      // src.convertTo(src, CV_8UC1); //直接转换会出错
      cv::cvtColor(src, src, cv::COLOR_GRAY2BGR); // 转换为彩色图
      // cv::imwrite(fileName + "src.png", src);
      // 步骤2：图像预处理
      int64_t begin =
      std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now().time_since_epoch())
      .count();
      cv::Mat gray, blurred, edges;
      cv::cvtColor(src, gray, cv::COLOR_BGR2GRAY); // 转换为灰度图
      PreProcessing(gray, blurred);
      /*边界检测分析：目标L将圆球内部及非目标物体内的边界去除
      sobel 噪声边界更少
      1 0 -1
      2 0 -2
      1 0 -1
      Scharr
      -3 0 3
      -10 0 10
      -3 0 3
      Laplace函数实现的方法是先用Sobel 算子计算二阶x和y导数，再求和
      */
      // cv::Mat dx, dy;
      // Sobel(gray, dx, CV_16S, 1, 0, 3, 1, 0, cv::BORDER_REPLICATE);
      // Sobel(gray, dy, CV_16S, 0, 1, 3, 1, 0, cv::BORDER_REPLICATE);
      // Canny(dx, dy, edges, 25, 50, false); // 边缘检测:高于200==strong edge,直接保留; 低于100，直接剔除; 100～
      // cv::imwrite("Compare/" + fileName + "graydx.png", dx);
      // cv::imwrite("Compare/" + fileName + "graydy.png", dy);
      // cv::imwrite("Compare/" + fileName + "grayedges.png", edges);
      // // Scharr(gray, dx, CV_16S, 1, 0);
      // // Scharr(gray, dy, CV_16S, 0, 1);
      // Laplacian(gray, dx, CV_16S, 3);
      // Laplacian(gray, dy, CV_16S, 3);
      // Canny(dx, dy, edges, 25, 50, false);
      // cv::imwrite("Compare/" + fileName + "blurreddx.png", dx);
      // cv::imwrite("Compare/" + fileName + "blurreddy.png", dy);
      // cv::imwrite("Compare/" + fileName + "blurrededges.png", edges);
      // 步骤3：寻找轮廓
      // 使用Hough变换检测圆
      /*
      src: 输入图像，必须是灰度图（8位）
      circles: 检测到的圆输出向量。每个圆由3个浮点数表示：(x, y, r)，其中(x,
      y)是圆心的坐标，r是圆的半径。
      cv::HOUGH_GRADIENT:
      这是用于圆检测的方法。目前OpenCV只实现了这种方法，它是一种利用边缘梯度进行圆心和半径检测的霍夫变换方法。
      1: 这是dp，累加器分辨率与图像分辨率的比率。在这个例子中，累加器的分辨率与输入图像相同。
      10:
      这是minDist参数，表示检测到的圆心之间的最小距离。在这里，它被设置为图像行数的八分之一，意味着检测到的圆的中心必须至少相距图像高度的八分之一，这样可以避免检测到多个相邻的圆。
      50:
      这是param1参数，用于Canny边缘检测器的高阈值。较低的阈值是这个值的一半，即100。这影响了圆检测的边缘阶段，边缘检测越敏感，可以检测到的圆就越多，但也可能导致更多的假阳性。
      cv::Canny(blurred, edges, 20, 50); // 边缘检测:高于200==strong edge,直接保留; 低于100，直接剔除; 100～200之间==weak edge：与strong edge则保留
      即此时小于param1/2.0的像素被剔除边界选择
      30:
      这是param2参数，累加器阈值。只有那些累加器中的值高于这个阈值的圆才会被考虑。较小的值意味着更多的圆会被检测到，但也可能包括一些假圆。
      累加器中的值 = 像素点个数:
      对于每一个边缘点 (x, y)，在累加器空间中对所有可能的圆心 (a, b) 和半径 r 进行投票。具体来说，对于每一个可能的半径 r，计算可能的圆心 (a, b) 并在累加器中相应的位置加一
      3: 这是minRadius，圆半径的最小值。在这个例子中，它被设置为0，意味着不对圆的最小半径设限。
      30: 这是maxRadius，圆半径的最大值。同样地，它被设置为0，意味着不对圆的最大半径设限。
      @param3:==3
      使用sobel进行水平及竖直梯度计算的kernelSize;代码默认参数，无法设置
      1 0 -1
      2 0 -2
      1 0 -1
      */
      std::vector<cv::Vec3f> circles;
      // cv::HoughCircles(blurred, circles, cv::HOUGH_GRADIENT, 1, 10, 50, 30, 5, 40);
      MyHoughCircles(blurred, circles, 1.0, 10, 50, 30, 5, 40);
      // 绘制检测圆
      for (size_t i = 0; i < circles.size(); i++)
      {
      cv::circle(src, cv::Point(circles[i][0], circles[i][1]), circles[i][2], cv::Scalar(255, 0, 0), 3); // 绘制圆
      }
      // std::cout << fileName << ": number of circles before deleted: " << circles.size() << std::endl;
      // cv::imwrite(fileName + "HoughCircles.png", src);
      int sumNumber = 100;
      if (circles.size() > sumNumber)
      {
      std::vector<CircleCenter> allDistance;
      int distance = 0;
      for (int i = 0; i < circles.size(); i++)
      {
      distance = 0;
      for (int j = 0; j < circles.size(); j++)
      {
      // 不加根号会放大远距离点的比重，导致某些点到最远点距离相近的点都被识别为目标点
      distance += sqrt((pow(circles[i][0] - circles[j][0], 2) + pow(circles[i][1] - circles[j][1], 2)));
      }
      CircleCenter temp;
      temp.circle = circles[i];
      temp.distance = distance;
      allDistance.push_back(temp);
      // std::cout << cvRound(circles[i][0]) << " " << cvRound(circles[i][1]) << " " << temp.distance << std::endl;
      }
      std::sort(allDistance.begin(), allDistance.end(), mycompare);
      circles.resize(sumNumber);
      for (int i = 0; i < sumNumber; i++)
      {
      circles[i] = allDistance[i].circle;
      }
      // 绘制被剔除圆
      for (size_t i = sumNumber; i < allDistance.size(); i++)
      {
      cv::circle(src, cv::Point(allDistance[i].circle[0], allDistance[i].circle[1]), allDistance[i].circle[2], cv::Scalar(0, 0, 255), 3);
      }
      }
      int64_t end =
      std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now().time_since_epoch())
      .count();
      std::cout << "Time Used:" << (double)(end - begin) / 1000 << "ms\n";
      std::cout << fileName << ": number of circles after deleted: " << circles.size() << std::endl;
      cv::imshow("Detected Ellipses", src);
      cv::imwrite("Result/" + fileName + "ellipse.png", src);
      return 0;
      }
      int main(int argc, char **argv)
      {
      QString fileType = ".dcm";
      QString dirpath(argv[1]);
      qDebug() << dirpath;
      QStringList roms_files;
      QDir dirtmp(dirpath);
      // 设置文件过滤器
      QStringList nameFilters;
      // 设置文件过滤格式
      nameFilters << "*" + fileType;
      // 将过滤后的文件名称存入到files列表中
      roms_files = dirtmp.entryList(nameFilters, QDir::Files | QDir::Readable, QDir::Name);
      qDebug() << roms_files;
      std::vector<std::string> vecFiles;
      if (roms_files.size() > 0)
      {
      vecFiles.reserve(static_cast<size_t>(roms_files.size()));
      for (int i = 0; i < roms_files.size(); i++)
      {
      int nPos = roms_files[i].indexOf(fileType);
      QString str = roms_files[i].mid(0, nPos);
      vecFiles.push_back(str.toStdString());
      }
      }
      for (int i = 0; i < vecFiles.size(); i++)
      {
      HoughCirclesExample(vecFiles[i], dirpath.toStdString(), fileType.toStdString());
      }
      // HoughCirclesExample("1", "0614");
      cvWaitKey(0);
      }
      \ No newline at end of file
      myHough.h 0 → 100644
      View file @ 3676ea7d
      #include <opencv2/opencv.hpp>
      using namespace cv;
      class MyNZPointSet
      {
      private:
      MyNZPointSet(const MyNZPointSet &other); // non-copyable
      public:
      Mat_<uchar> positions;
      MyNZPointSet(int rows, int cols) : positions(rows, cols, (uchar)0)
      {
      }
      void insert(const Point &pt)
      {
      positions(pt) = 1;
      }
      void insert(const MyNZPointSet &from)
      {
      cv::bitwise_or(from.positions, positions, positions);
      }
      };
      class MyHoughCirclesAccumInvoker : public ParallelLoopBody
      {
      public:
      MyHoughCirclesAccumInvoker(const Mat &_edges, const Mat &_dx, const Mat &_dy, int _minRadius, int _maxRadius, float _idp,
      std::vector<Mat> &_accumVec, MyNZPointSet &_nz, Mutex &_mtx) : edges(_edges), dx(_dx), dy(_dy), minRadius(_minRadius), maxRadius(_maxRadius), idp(_idp),
      accumVec(_accumVec), nz(_nz), mutex(_mtx)
      {
      acols = cvCeil(edges.cols * idp), arows = cvCeil(edges.rows * idp);
      astep = acols + 2;
      }
      ~MyHoughCirclesAccumInvoker() {}
      void operator()(const Range &boundaries) const CV_OVERRIDE
      {
      Mat accumLocal = Mat(arows + 2, acols + 2, CV_32SC1, Scalar::all(0));
      int *adataLocal = accumLocal.ptr<int>();
      MyNZPointSet nzLocal(nz.positions.rows, nz.positions.cols);
      int startRow = boundaries.start;
      int endRow = boundaries.end;
      int numCols = edges.cols;
      if (edges.isContinuous() && dx.isContinuous() && dy.isContinuous())
      {
      numCols *= (boundaries.end - boundaries.start);
      endRow = boundaries.start + 1;
      }
      // Accumulate circle evidence for each edge pixel
      for (int y = startRow; y < endRow; ++y)
      {
      const uchar *edgeData = edges.ptr<const uchar>(y);
      const short *dxData = dx.ptr<const short>(y);
      const short *dyData = dy.ptr<const short>(y);
      int x = 0;
      for (; x < numCols; ++x)
      {
      #if CV_SIMD
      {
      v_uint8 v_zero = vx_setzero_u8();
      for (; x <= numCols - 2 * v_uint8::nlanes; x += 2 * v_uint8::nlanes)
      {
      v_uint8 v_edge1 = (vx_load(edgeData + x) != v_zero);
      v_uint8 v_edge2 = (vx_load(edgeData + x + v_uint8::nlanes) != v_zero);
      if (v_check_any(v_edge1))
      {
      x += v_scan_forward(v_edge1);
      goto _next_step;
      }
      if (v_check_any(v_edge2))
      {
      x += v_uint8::nlanes + v_scan_forward(v_edge2);
      goto _next_step;
      }
      }
      }
      #endif
      for (; x < numCols && !edgeData[x]; ++x)
      ;
      if (x == numCols)
      continue;
      #if CV_SIMD
      _next_step:
      #endif
      float vx, vy;
      int sx, sy, x0, y0, x1, y1;
      vx = dxData[x];
      vy = dyData[x];
      if (vx == 0 && vy == 0)
      continue;
      float mag = std::sqrt(vx * vx + vy * vy);
      if (mag < 1.0f)
      continue;
      Point pt = Point(x % edges.cols, y + x / edges.cols);
      nzLocal.insert(pt);
      sx = cvRound((vx * idp) * 1024 / mag);
      sy = cvRound((vy * idp) * 1024 / mag);
      x0 = cvRound((pt.x * idp) * 1024);
      y0 = cvRound((pt.y * idp) * 1024);
      // Step from min_radius to max_radius in both directions of the gradient
      for (int k1 = 0; k1 < 2; k1++)
      {
      x1 = x0 + minRadius * sx;
      y1 = y0 + minRadius * sy;
      for (int r = minRadius; r <= maxRadius; x1 += sx, y1 += sy, r++)
      {
      int x2 = x1 >> 10, y2 = y1 >> 10;
      if ((unsigned)x2 >= (unsigned)acols ||
      (unsigned)y2 >= (unsigned)arows)
      break;
      adataLocal[y2 * astep + x2]++;
      }
      sx = -sx;
      sy = -sy;
      }
      }
      }
      { // TODO Try using TLSContainers
      AutoLock lock(mutex);
      accumVec.push_back(accumLocal);
      nz.insert(nzLocal);
      }
      }
      private:
      const Mat &edges, &dx, &dy;
      int minRadius, maxRadius;
      float idp;
      std::vector<Mat> &accumVec;
      MyNZPointSet &nz;
      int acols, arows, astep;
      Mutex &mutex;
      };
      // 查找4领域内的最大值同时符合阈值的累加器结果
      class MyHoughCirclesFindCentersInvoker : public ParallelLoopBody
      {
      public:
      MyHoughCirclesFindCentersInvoker(const Mat &_accum, std::vector<int> &_centers, int _accThreshold, Mutex &_mutex) : accum(_accum), centers(_centers), accThreshold(_accThreshold), _lock(_mutex)
      {
      acols = accum.cols;
      arows = accum.rows;
      adata = accum.ptr<int>();
      }
      ~MyHoughCirclesFindCentersInvoker() {}
      void operator()(const Range &boundaries) const CV_OVERRIDE
      {
      int startRow = boundaries.start;
      int endRow = boundaries.end;
      std::vector<int> centersLocal;
      bool singleThread = (boundaries == Range(1, accum.rows - 1));
      startRow = max(1, startRow);
      endRow = min(arows - 1, endRow);
      // Find possible circle centers
      for (int y = startRow; y < endRow; ++y)
      {
      int x = 1;
      int base = y * acols + x;
      for (; x < acols - 1; ++x, ++base)
      {
      if (adata[base] > accThreshold &&
      adata[base] > adata[base - 1] && adata[base] >= adata[base + 1] &&
      adata[base] > adata[base - acols] && adata[base] >= adata[base + acols])
      centersLocal.push_back(base);
      }
      }
      if (!centersLocal.empty())
      {
      if (singleThread)
      centers = centersLocal;
      else
      {
      AutoLock alock(_lock);
      centers.insert(centers.end(), centersLocal.begin(), centersLocal.end());
      }
      }
      }
      private:
      const Mat &accum;
      std::vector<int> &centers;
      int accThreshold;
      int acols, arows;
      const int *adata;
      Mutex &_lock;
      };
      struct EstimatedCircle
      {
      EstimatedCircle(Vec3f _c, int _accum) : c(_c), accum(_accum) {}
      Vec3f c;
      int accum;
      };
      struct hough_cmp_gt
      {
      hough_cmp_gt(const int *_aux) : aux(_aux) {}
      inline bool operator()(int l1, int l2) const
      {
      return aux[l1] > aux[l2] || (aux[l1] == aux[l2] && l1 < l2);
      }
      const int *aux;
      };
      static bool cmpAccum(const EstimatedCircle &left, const EstimatedCircle &right)
      {
      // Compare everything so the order is completely deterministic
      // Larger accum first
      if (left.accum > right.accum)
      return true;
      else if (left.accum < right.accum)
      return false;
      // Larger radius first
      else if (left.c[2] > right.c[2])
      return true;
      else if (left.c[2] < right.c[2])
      return false;
      // Smaller X
      else if (left.c[0] < right.c[0])
      return true;
      else if (left.c[0] > right.c[0])
      return false;
      // Smaller Y
      else if (left.c[1] < right.c[1])
      return true;
      else if (left.c[1] > right.c[1])
      return false;
      // Identical - neither object is less than the other
      else
      return false;
      }
      template <typename T>
      static bool CheckDistance(const std::vector<T> &circles, size_t endIdx, const T &circle, float minDist2)
      {
      bool goodPoint = true;
      for (uint j = 0; j < endIdx; ++j)
      {
      T pt = circles[j];
      float distX = circle[0] - pt[0], distY = circle[1] - pt[1];
      if (distX * distX + distY * distY < minDist2)
      {
      goodPoint = false;
      break;
      }
      }
      return goodPoint;
      }
      template <typename T>
      static void RemoveOverlaps(std::vector<T> &circles, float minDist)
      {
      if (circles.size() <= 1u)
      return;
      float minDist2 = minDist * minDist;
      size_t endIdx = 1;
      for (size_t i = 1; i < circles.size(); ++i)
      {
      T circle = circles[i];
      if (CheckDistance(circles, endIdx, circle, minDist2))
      {
      circles[endIdx] = circle;
      ++endIdx;
      }
      }
      circles.resize(endIdx);
      }
      static inline Vec3f GetCircle(const EstimatedCircle &est)
      {
      return est.c;
      }
      static void CreateCircles(const std::vector<EstimatedCircle> &circlesEst, std::vector<Vec3f> &circles)
      {
      std::transform(circlesEst.begin(), circlesEst.end(), std::back_inserter(circles), GetCircle);
      }
      template <class NZPoints>
      class MyHoughCircleEstimateRadiusInvoker : public ParallelLoopBody
      {
      public:
      MyHoughCircleEstimateRadiusInvoker(const NZPoints &_nz, int _nzSz, const std::vector<int> &_centers, std::vector<EstimatedCircle> &_circlesEst,
      int _acols, int _accThreshold, int _minRadius, int _maxRadius,
      float _dp, Mutex &_mutex) : nz(_nz), nzSz(_nzSz), centers(_centers), circlesEst(_circlesEst), acols(_acols), accThreshold(_accThreshold),
      minRadius(_minRadius), maxRadius(_maxRadius), dr(_dp), _lock(_mutex)
      {
      minRadius2 = (float)minRadius * minRadius;
      maxRadius2 = (float)maxRadius * maxRadius;
      centerSz = (int)centers.size();
      CV_Assert(nzSz > 0);
      }
      ~MyHoughCircleEstimateRadiusInvoker() {}
      protected:
      inline int filterCircles(const Point2f &curCenter, float *ddata) const;
      void operator()(const Range &boundaries) const CV_OVERRIDE
      {
      std::vector<EstimatedCircle> circlesLocal;
      const int nBinsPerDr = 10; // 于该参数无关
      // const int nBinsPerDr = 1;
      int nBins = cvRound((maxRadius - minRadius) / dr * nBinsPerDr);
      AutoBuffer<int> bins(nBins);
      AutoBuffer<float> distBuf(nzSz), distSqrtBuf(nzSz);
      float *ddata = distBuf.data();
      float *dSqrtData = distSqrtBuf.data();
      bool singleThread = (boundaries == Range(0, centerSz));
      int i = boundaries.start;
      // For each found possible center
      // Estimate radius and check support
      for (; i < boundaries.end; ++i)
      {
      int ofs = centers[i];
      int y = ofs / acols;
      int x = ofs - y * acols;
      // Calculate circle's center in pixels
      Point2f curCenter = Point2f((x + 0.5f) * dr, (y + 0.5f) * dr);
      int nzCount = filterCircles(curCenter, ddata);
      int maxCount = 0;
      float rBest = 0;
      if (nzCount)
      {
      Mat_<float> distMat(1, nzCount, ddata);
      Mat_<float> distSqrtMat(1, nzCount, dSqrtData);
      sqrt(distMat, distSqrtMat);
      memset(bins.data(), 0, sizeof(bins[0]) * bins.size());
      for (int k = 0; k < nzCount; k++)
      {
      int bin = std::max(0, std::min(nBins - 1, cvRound((dSqrtData[k] - minRadius) / dr * nBinsPerDr)));
      bins[bin]++;
      }
      for (int j = nBins - 1; j > 0; j--)
      {
      if (bins[j])
      {
      int upbin = j;
      int curCount = 0;
      for (; j > upbin - nBinsPerDr && j >= 0; j--)
      {
      curCount += bins[j];
      }
      float rCur = (upbin + j) / 2.f / nBinsPerDr * dr + minRadius;
      if ((curCount * rBest >= maxCount * rCur) ||
      (rBest < FLT_EPSILON && curCount >= maxCount))
      {
      rBest = rCur;
      maxCount = curCount;
      }
      }
      }
      }
      // Check if the circle has enough support
      if (maxCount > accThreshold)
      {
      circlesLocal.push_back(EstimatedCircle(Vec3f(curCenter.x, curCenter.y, rBest), maxCount));
      }
      }
      if (!circlesLocal.empty())
      {
      std::sort(circlesLocal.begin(), circlesLocal.end(), cmpAccum);
      if (singleThread)
      {
      std::swap(circlesEst, circlesLocal);
      }
      else
      {
      AutoLock alock(_lock);
      if (circlesEst.empty())
      std::swap(circlesEst, circlesLocal);
      else
      circlesEst.insert(circlesEst.end(), circlesLocal.begin(), circlesLocal.end());
      }
      }
      }
      private:
      const NZPoints &nz;
      int nzSz;
      const std::vector<int> &centers;
      std::vector<EstimatedCircle> &circlesEst;
      int acols, accThreshold, minRadius, maxRadius;
      float dr;
      int centerSz;
      float minRadius2, maxRadius2;
      Mutex &_lock;
      };
      // 在curCenter的最大查找半径区域内遍历点，将点到curCenter的距离满足最大最小半径的点存入ddata
      template <>
      inline int MyHoughCircleEstimateRadiusInvoker<MyNZPointSet>::filterCircles(const Point2f &curCenter, float *ddata) const
      {
      int nzCount = 0;
      const Mat_<uchar> &positions = nz.positions;
      const int rOuter = maxRadius + 1;
      const Range xOuter = Range(std::max(int(curCenter.x - rOuter), 0), std::min(int(curCenter.x + rOuter), positions.cols));
      const Range yOuter = Range(std::max(int(curCenter.y - rOuter), 0), std::min(int(curCenter.y + rOuter), positions.rows));
      #if CV_SIMD
      float v_seq[v_float32::nlanes];
      for (int i = 0; i < v_float32::nlanes; ++i)
      v_seq[i] = (float)i;
      const v_float32 v_minRadius2 = vx_setall_f32(minRadius2);
      const v_float32 v_maxRadius2 = vx_setall_f32(maxRadius2);
      const v_float32 v_curCenterX_0123 = vx_setall_f32(curCenter.x) - vx_load(v_seq);
      #endif
      for (int y = yOuter.start; y < yOuter.end; y++)
      {
      const uchar *ptr = positions.ptr(y, 0);
      float dy = curCenter.y - y;
      float dy2 = dy * dy;
      int x = xOuter.start;
      #if CV_SIMD
      {
      const v_float32 v_dy2 = vx_setall_f32(dy2);
      const v_uint32 v_zero_u32 = vx_setall_u32(0);
      float CV_DECL_ALIGNED(CV_SIMD_WIDTH) rbuf[v_float32::nlanes];
      int CV_DECL_ALIGNED(CV_SIMD_WIDTH) rmask[v_int32::nlanes];
      for (; x <= xOuter.end - v_float32::nlanes; x += v_float32::nlanes)
      {
      v_uint32 v_mask = vx_load_expand_q(ptr + x);
      v_mask = v_mask != v_zero_u32;
      v_float32 v_x = v_cvt_f32(vx_setall_s32(x));
      v_float32 v_dx = v_x - v_curCenterX_0123;
      v_float32 v_r2 = (v_dx * v_dx) + v_dy2;
      v_float32 vmask = (v_minRadius2 <= v_r2) & (v_r2 <= v_maxRadius2) & v_reinterpret_as_f32(v_mask);
      if (v_check_any(vmask))
      {
      v_store_aligned(rmask, v_reinterpret_as_s32(vmask));
      v_store_aligned(rbuf, v_r2);
      for (int i = 0; i < v_int32::nlanes; ++i)
      if (rmask[i])
      ddata[nzCount++] = rbuf[i];
      }
      }
      }
      #endif
      for (; x < xOuter.end; x++)
      {
      if (ptr[x])
      {
      float _dx = curCenter.x - x;
      float _r2 = _dx * _dx + dy2;
      if (minRadius2 <= _r2 && _r2 <= maxRadius2)
      {
      ddata[nzCount++] = _r2;
      }
      }
      }
      }
      return nzCount;
      }
      void MyHoughCircles(cv::Mat src, std::vector<Vec3f> &circles, double dp, int minDist, int param1, int accThreshold, int minRadius, int maxRadius)
      {
      cv::Mat accumSrc, centerSrc;
      cv::cvtColor(src, accumSrc, cv::COLOR_GRAY2BGR);
      cv::cvtColor(src, centerSrc, cv::COLOR_GRAY2BGR);
      cv::Mat edges, dx, dy;
      Sobel(src, dx, CV_16S, 1, 0, 3, 1, 0, cv::BORDER_REPLICATE);
      Sobel(src, dy, CV_16S, 0, 1, 3, 1, 0, cv::BORDER_REPLICATE);
      Canny(dx, dy, edges, param1 / 2, param1, false);
      // 1.计算累加器结果
      cv::Mutex mtx;
      int numThreads = std::max(1, cv::getNumThreads());
      std::vector<cv::Mat> accumVec;
      MyNZPointSet nz(src.rows, src.cols);
      cv::parallel_for_(Range(0, edges.rows),
      MyHoughCirclesAccumInvoker(edges, dx, dy, minRadius, maxRadius, 1.0, accumVec, nz, mtx),
      numThreads);
      int nzSz = cv::countNonZero(nz.positions);
      cv::Mat accum = accumVec[0];
      for (size_t i = 1; i < accumVec.size(); i++)
      {
      accum += accumVec[i];
      }
      // 累加器结果
      for (int i = 0; i < accum.rows; i++)
      {
      for (int j = 0; j < accum.cols; j++)
      {
      if (accum.at<int>(i, j) > accThreshold)
      // if (1)
      {
      // std::cout << i << " , " << j << " : " << accum.at<int>(i, j) << std::endl;
      // showSrc.at<cv::Vec3b>(i, j) = cv::Vec3b(accum.at<int>(i, j), 0, 0); // blue green red
      accumSrc.at<cv::Vec3b>(i, j) = cv::Vec3b(255, 0, 0); // blue green red
      }
      }
      }
      // cv::imwrite("accumSrc.png", accumSrc);
      accumVec.clear();
      // 2.分析累加器中的潜在center
      std::vector<int> centers;
      // 4 rows when multithreaded because there is a bit overhead
      // and on the other side there are some row ranges where centers are concentrated
      cv::parallel_for_(Range(1, accum.rows - 1),
      MyHoughCirclesFindCentersInvoker(accum, centers, accThreshold, mtx),
      (numThreads > 1) ? ((accum.rows - 2) / 4) : 1);
      int centerCnt = (int)centers.size();
      // std::cout << "centerCnt " << centerCnt << std::endl;
      for (int i = 0; i < centerCnt; i++)
      {
      accumSrc.at<cv::Vec3b>(centers[i] / accum.cols, centers[i] % accum.cols) = cv::Vec3b(0, 0, 255); // 访问为先行号(Y轴)，后列号(X轴)
      }
      // cv::imwrite("centerSrc.png", accumSrc);
      // 3.分析潜在center的实际center与radius
      if (centerCnt == 0)
      return;
      std::sort(centers.begin(), centers.end(), hough_cmp_gt(accum.ptr<int>()));
      // 只有28个点符合要求
      std::vector<EstimatedCircle> circlesEst;
      parallel_for_(Range(0, centerCnt),
      MyHoughCircleEstimateRadiusInvoker<MyNZPointSet>(nz, nzSz, centers, circlesEst, accum.cols,
      accThreshold / 2, minRadius, maxRadius, dp, mtx),
      numThreads);
      // std::cout << "circlesEst " << circlesEst.size() << std::endl;
      for (int i = 0; i < circlesEst.size(); i++)
      {
      centerSrc.at<cv::Vec3b>(circlesEst[i].c[1], circlesEst[i].c[0]) = cv::Vec3b(0, 0, circlesEst[i].accum); // 访问为先行号(Y轴)，后列号(X轴)
      }
      std::sort(circlesEst.begin(), circlesEst.end(), cmpAccum);
      // Create Circles
      CreateCircles(circlesEst, circles);
      RemoveOverlaps(circles, minDist);
      // std::cout << "RemoveOverlaps " << circles.size() << std::endl;
      // cv::imwrite("fiannlCenterSrc.png", centerSrc);
      }
      \ No newline at end of file
      Markdown is supported
      0% or
      You are about to add 0 people to the discussion. Proceed with caution.
      Finish editing this message first!
      Please register or to comment