title: Detect clusters of circular objects by iterative adaptive thresholding and shape analysis
body: <p>I have been developing an application to count circular objects such as bacterial colonies from pictures. </p>
  
  <p>What make it easy is the fact that the objects are generally well distinct from the background.</p>
  
  <p>However, few difficulties make the analysis tricky:</p>
  
  <ol>
  <li>The background will present gradual as well as rapid intensity change.</li>
  <li>In the edges of the container, the object will be elliptic rather than circular.</li>
  <li>The edges of the objects are sometimes rather fuzzy.</li>
  <li>The objects will cluster.</li>
  <li>The object can be very small (6px of diameter)</li>
  <li>Ultimately, the algorithms will be used (via GUI) by people that do not have deep understanding of image analysis, so the parameters must be intuitive and very few.</li>
  </ol>
  
  <p>The problem has been address many times in the scientific literature and "solved", for instance, using circular Hough transform or watershed approaches, but I have never been satisfied by the results.</p>
  
  <p>One simple approach that was described is to get the foreground by adaptive thresholding and split (as I described in <a href="https://stackoverflow.com/questions/10313602/reshaping-noisy-coin-into-a-circle-form/10324759#10324759">this post</a>) the clustered objects using distance transform.</p>
  
  <p>I have successfully implemented this method, but it could not always deal with sudden change in intensity. Also, I have been asked by peers to come out with a more "novel" approach.</p>
  
  <p><strong>I therefore was looking for a new method to extract foreground.</strong></p>
  
  <p>I therefore investigated other thresholding/blob detection methods.
  I tried MSERs but found out that they were not very robust and quite slow in my case.</p>
  
  <p>I eventually came out with an algorithm that, so far, gives me excellent results:</p>
  
  <ol>
  <li>I split the three channels of my image and reduce their noise (blur/median blur). For each channel:</li>
  <li>I apply a manual implementation of the first step of adaptive thresholding by calculating the absolute difference between the original channel and a convolved (by a large kernel blur) one. Then, for all the relevant values of threshold:</li>
  <li>I apply a threshold on the result of 2)</li>
  <li>find contours</li>
  <li>validate or invalidate contours on the grant of their shape (size, area, convexity...)</li>
  <li>only the valid continuous regions (<em>i.e.</em> delimited by contours) are then redrawn in an accumulator (1 accumulator per channel).</li>
  <li>After accumulating continuous regions over values of threshold, I end-up with a map of "scores of regions". The regions with the highest intensity being those that fulfilled the the morphology filter criteria the most often.</li>
  <li>The three maps (one per channel) are then converted to grey-scale and thresholded (the threshold is controlled by the user)</li>
  </ol>
  
  <p>Just to show you the kind of image I have to work with: 
  <img src="https://i.stack.imgur.com/vXjUt.jpg" alt="enter image description here">
  This picture represents part of 3 sample images in the top and the result of my algorithm (blue = foreground) of the respective parts in the bottom.</p>
  
  <p>Here is my C++ implementation of : 3-7</p>
  
  <pre><code>/*
   * cv::Mat dst[3] is the result of the absolute difference between original and convolved channel.
   * MCF(std::vector&lt;cv::Point&gt;, int, int) is a filter function that returns an positive int only if the input contour is valid.
   */
  
  /* Allocate 3 matrices (1 per channel)*/
  cv::Mat accu[3];
  
  /* We define the maximal threshold to be tried as half of the absolute maximal value in each channel*/
  int maxBGR[3];
  for(unsigned int i=0; i&lt;3;i++){
      double min, max;
      cv::minMaxLoc(dst[i],&amp;min,&amp;max);
      maxBGR[i] = max/2;
      /* In addition, we fill accumulators by zeros*/
      accu[i]=cv::Mat(compos[0].rows,compos[0].cols,CV_8U,cv::Scalar(0));
  }
  /* This loops are intended to be multithreaded using
  #pragma omp parallel for collapse(2) schedule(dynamic)
  For each channel */
  for(unsigned int i=0; i&lt;3;i++){
      /* For each value of threshold (m_step can be &gt; 1 in order to save time)*/
      for(int j=0;j&lt;maxBGR[i] ;j += m_step ){
              /* Temporary matrix*/
              cv::Mat tmp;
              std::vector&lt;std::vector&lt;cv::Point&gt; &gt; contours;
              /* Thresholds dst by j*/
              cv::threshold(dst[i],tmp, j, 255, cv::THRESH_BINARY);
              /* Finds continous regions*/
              cv::findContours(tmp, contours, CV_RETR_LIST, CV_CHAIN_APPROX_TC89_L1);
              if(contours.size() &gt; 0){
                  /* Tests each contours*/
                  for(unsigned int k=0;k&lt;contours.size();k++){
                      int valid = MCF(contours[k],m_minRad,m_maxRad);
                      if(valid&gt;0){
                          /* I found that redrawing was very much faster if the given contour was copied in a smaller container.
                           * I do not really understand why though. For instance,
                           cv::drawContours(miniTmp,contours,k,cv::Scalar(1),-1,8,cv::noArray(), INT_MAX, cv::Point(-rect.x,-rect.y));
                           is slower especially if contours is very long.
                           */
                          std::vector&lt;std::vector&lt;cv::Point&gt; &gt; tpv(1);
                          std::copy(contours.begin()+k, contours.begin()+k+1, tpv.begin());
                          /* We make a Roi here*/
                          cv::Rect rect = cv::boundingRect(tpv[0]);
                          cv::Mat miniTmp(rect.height,rect.width,CV_8U,cv::Scalar(0));
                          cv::drawContours(miniTmp,tpv,0,cv::Scalar(1),-1,8,cv::noArray(), INT_MAX, cv::Point(-rect.x,-rect.y));
                          accu[i](rect)    = miniTmp + accu[i](rect);
                      }
                  }
              }
          }
      }
  /* Make the global scoreMap*/
  cv::merge(accu,3,scoreMap);
  /* Conditional noise removal*/
  if(m_minRad&gt;2)
      cv::medianBlur(scoreMap,scoreMap,3);
  cvtColor(scoreMap,scoreMap,CV_BGR2GRAY);
  </code></pre>
  
  <p>I have two questions:</p>
  
  <ol>
  <li><p>What is the name of such foreground extraction approach and do you see any reason for which it could be improper to use it in this case ?</p></li>
  <li><p>Since recursively finding and drawing contours is quite intensive, I would like to make my algorithm faster. Can you indicate me any way to achieve this goal ?</p></li>
  </ol>
  
  <p>Thank you very much for you help,</p>