18

I want to be able to find points in images that are the centre of a radial gradient like the one shown in the left picture below. Any ideas on how I could use a Hough transform or some other computer vision method?

Thanks

enter image description here

example search image:

enter image description here

Royi
  • 33,983
  • 4
  • 72
  • 179
waspinator
  • 743
  • 1
  • 7
  • 12
  • Great question! – Spacey Aug 01 '12 at 18:11
  • Also, take a look at Roberts' Cross: (http://en.wikipedia.org/wiki/Roberts_Cross) as an example of a way to estimate gradients. – Spacey Aug 01 '12 at 19:25
  • looks like a smaller sobel operator. I'm not sure how to use that to find a radial gradient though – waspinator Aug 01 '12 at 19:32
  • @waspinator: well have you run a sobel operator on your image and looked at the output? It's like the 2D equivalent of taking the derivative of a 1D function, so it should cross 0 at the local minima or maxima? – endolith Aug 01 '12 at 19:33
  • If you want to find a point on the top/bottom of a peak/valley, then you can use an iterative algorithm like gradient ascent/descent. The image in your case is actually the cost function. At every point a 2-D gradient ascent/descent will calculate the gradient and move along it accordingly. – Spacey Aug 01 '12 at 19:47
  • output of a sobel operator using matlab http://i.imgur.com/gP1gd.png – waspinator Aug 01 '12 at 19:49
  • gradient ascent/descent would only find a single local maximum, and would also find gradients which aren't radial, which I wouldn't want. example search image: http://i.imgur.com/ToxhL.png – waspinator Aug 01 '12 at 19:59
  • 1
    For a simple Hough-like approach that would probably work you could try this: for every pixel of the image, calculate the gradient direction and render a short line segment in the direction of the gradient starting at this pixel into an accumulator. The center points you are looking for should be the highest peaks in the accumulator (by a large margin). – koletenbert Aug 01 '12 at 21:18
  • @waspinator Ok, perhaps you can add that image to the question so that we have a clearer scope. Thanks! – Spacey Aug 01 '12 at 21:57
  • Hi, I know it's been a while but I ran into the exact same problem. i.e. finding circular gradients in an image. Did any of the answers below solve your issue, or did you find another way? – Cape Code Jan 14 '14 at 20:16

3 Answers3

7

I was working in opencv and trying to find the peak of a gradient generated by distance transform. I realised that using morphologic operations (erosion/dilatation) in grey-scsale images was very useful in this case. If you erode dilate a grey-scale image, any pixel is going to take the value of the lower/highest neighbour. You can therefore find peaks of intensity in gradients by subtracting grey-scale image from the same dilated/eroded image. Here is my result: enter image description here

And a way to do it in OpenCV/Cpp:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

int main( int argc, char** argv ){

    cv::Mat objects, img ,peaks,BGR;
    std::vector<std::vector<cv::Point> > contours;
    /* Reads the image*/
    BGR=cv::imread(argv[1]);
    /* Converts it to Grayscale*/
    cv::cvtColor(BGR,img,CV_BGR2GRAY);
    /* Devine where are the objects*/
    cv::threshold(img,objects,0,255,cv::THRESH_BINARY);
    /* In order to find the local maxima, "distance"
     * is subtracted from the result of the dilatation of
     * "distance". All the peaks keep the save value */
    cv::dilate(img,peaks,cv::Mat(),cv::Point(-1,-1),3);
    cv::dilate(objects,objects,cv::Mat(),cv::Point(-1,-1),3);

    /* Now all the peaks should be exactely 0*/
    peaks=peaks-img;

    /* And the non-peaks 255*/
    cv::threshold(peaks,peaks,0,255,cv::THRESH_BINARY);
    peaks.convertTo(peaks,CV_8U);

    /* Only the zero values of "peaks" that are non-zero
     * in "objects" are the real peaks*/
    cv::bitwise_xor(peaks,objects,peaks);

    /* The peaks that are distant from less than
     * 2 pixels are merged by dilatation */
    cv::dilate(peaks,peaks,cv::Mat(),cv::Point(-1,-1),1);

    /* In order to map the peaks, findContours() is used.
     * The results are stored in "contours" */
    cv::findContours(peaks, contours, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
    /* just draw them and save the image */
    cv::drawContours(BGR,contours,-1,cv::Scalar(255,0,0),-1);
    cv::imwrite("result.png",BGR);

    return 1;
}
5

Here is what I have so far. The way I'm populating my Hough space is far from optimal. I'm pretty sure there's some vectorization I can do to make it faster. I'm using Matlab R2011a. Original Image

Suggestions are appreciated, Thank you.

enter image description here

clear all; clc; close all;

%% read in image and find gradient information
img = rgb2gray(imread('123.png'));
[rows, columns] = size(img);
[dx, dy] = gradient(double(img));
[x y] = meshgrid(1:columns, 1:rows);
u = dx;
v = dy;
imshow(img);
hold on
quiver(x, y, u, v)


%% create Hough space and populate
hough_space = zeros(size(img));

for i = 1:columns
  for j = 1:rows

    X1 = i;
    Y1 = j;
    X2 = round(i + dx(j,i));
    Y2 = round(j + dy(j,i));
    increment = 1;

    slope = (Y2 - Y1) / (X2 - X1);
    y_intercept = Y1 - slope * X1;

    X3 = X1 + 5;

    if X3 < columns && X3 > 1
      Y3 = slope * X3 + y_intercept;
      if Y3 < rows && Y3 > 1
        hough_space = func_Drawline(hough_space, Y1, X1, floor(Y3), floor(X3), increment);
      end
    end
  end
end

imtool(hough_space)

I modified a draw line function I found on matlab central to increment by a pixel by a value instead of setting a pixel to a value

function Img = func_DrawLine(Img, X0, Y0, X1, Y1, nG)
% Connect two pixels in an image with the desired graylevel
%
% Command line
% ------------
% result = func_DrawLine(Img, X1, Y1, X2, Y2)
% input:    Img : the original image.
%           (X1, Y1), (X2, Y2) : points to connect.
%           nG : the gray level of the line.
% output:   result
%
% Note
% ----
%   Img can be anything
%   (X1, Y1), (X2, Y2) should be NOT be OUT of the Img
%
%   The computation cost of this program is around half as Cubas's [1]
%   [1] As for Cubas's code, please refer  
%   http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4177  
%
% Example
% -------
% result = func_DrawLine(zeros(5, 10), 2, 1, 5, 10, 1)
% result =
%      0     0     0     0     0     0     0     0     0     0
%      1     1     1     0     0     0     0     0     0     0
%      0     0     0     1     1     1     0     0     0     0
%      0     0     0     0     0     0     1     1     1     0
%      0     0     0     0     0     0     0     0     0     1
%
%
% Jing Tian Oct. 31 2000
% scuteejtian@hotmail.com
% This program is written in Oct.2000 during my postgraduate in 
% GuangZhou, P. R. China.
% Version 1.0

Img(X0, Y0) = Img(X0, Y0) + nG;
Img(X1, Y1) = Img(X1, Y1) + nG;
if abs(X1 - X0) <= abs(Y1 - Y0)
   if Y1 < Y0
      k = X1; X1 = X0; X0 = k;
      k = Y1; Y1 = Y0; Y0 = k;
   end
   if (X1 >= X0) & (Y1 >= Y0)
      dy = Y1-Y0; dx = X1-X0;
      p = 2*dx; n = 2*dy - 2*dx; tn = dy;
      while (Y0 < Y1)
         if tn >= 0
            tn = tn - p;
         else
            tn = tn + n; X0 = X0 + 1;
         end
         Y0 = Y0 + 1; Img(X0, Y0) = Img(X0, Y0) + nG;
      end
   else
      dy = Y1 - Y0; dx = X1 - X0;
      p = -2*dx; n = 2*dy + 2*dx; tn = dy;
      while (Y0 <= Y1)
         if tn >= 0
            tn = tn - p;
         else
            tn = tn + n; X0 = X0 - 1;
         end
         Y0 = Y0 + 1; Img(X0, Y0) = Img(X0, Y0) + nG;
      end
   end
else if X1 < X0
      k = X1; X1 = X0; X0 = k;
      k = Y1; Y1 = Y0; Y0 = k;
   end
   if (X1 >= X0) & (Y1 >= Y0)
      dy = Y1 - Y0; dx = X1 - X0;
      p = 2*dy; n = 2*dx-2*dy; tn = dx;
      while (X0 < X1)
         if tn >= 0
            tn = tn - p;
         else
            tn = tn + n; Y0 = Y0 + 1;
         end
         X0 = X0 + 1; Img(X0, Y0) = Img(X0, Y0) + nG;
      end
   else
      dy = Y1 - Y0; dx = X1 - X0;
      p = -2*dy; n = 2*dy + 2*dx; tn = dx;
      while (X0 < X1)
         if tn >= 0
            tn = tn - p;
         else
            tn = tn + n; Y0 = Y0 - 1;
         end
         X0 = X0 + 1; Img(X0, Y0) = Img(X0, Y0) + nG;
      end
   end
end
Rick Smith
  • 123
  • 5
waspinator
  • 743
  • 1
  • 7
  • 12
  • I think I will attribute the bounty to your answer, as nobody else bothered to contribute. It's not exactly what I want but it's the closest of the 3. Did you further improve this method? – Cape Code Jan 21 '14 at 20:46
1

Run a Histogram of Oriented Gradients over patches of the image - the peak in each of those histograms will give you the dominant direction of that patch (like the arrows you show).

Find where all those arrows intersect - if that point is inside the object it could be the centre of a radial gradient.

Martin Thompson
  • 1,432
  • 11
  • 17