Coordinates of every point on a circle's circumference

Question

First, please note, that this question is not a duplicate of these: 1st , 2nd , and 3rd.

I am using delphi and openCV, but I am looking for an algorithm, a solution regardless of the language.

For the purpose of a precise image analysis, I need to check for changes in pixel intensity in circular areas. So I read pixel values on a circumference of continuously growing circle. To be able to do that, I of course need to know coordinates of the pixels.

The best solution I found is y:= Round(centerY + radius * sin(angle)), x:= Round(centerX + radius * cos(angle)), while because counting with only 360 degrees is hardly enough, when the radius of the circle is larger, than circa 60px, the angle is being counted like this angle:= angle + (360 / (2 * 3.14 * currentRadius)) -> I sweep through every value from 0 to 360, while the value is being incremented by a fraction of 360/circumference of the circle in pixels. But this approach is not very precise. The larger the circle, the smaller the fraction of the angle needs to be and the precission suffers from the inaccuracy of Pi, plus the rounding.

If I use the mentioned method, and try to draw the counted pixels with this code:

  centerX:= 1700;
  centerY:= 1200;
  maxRadius:= 500;

  for currentRadius:= 80 to maxRadius do
  begin

    angle:= 0;
    while angle < 360 do
    begin

      xI:= Round(centerX + currentRadius * cos(angle));
      yI:= Round(centerY + currentRadius * sin(angle));
      angle:= angle + (360 / (2 * 3.14 * currentRadius));

      //this is openCV function, to test the code, you can use anything, that will draw a dot...
      cvLine(image,cvPoint(xI,yI),cvPoint(xI,yI),CV_RGB(0, 255, 0));

    end;

  end;

the result is this:

It is not bad, but taking into account, that rougly a third of all pixels in the circular area are black, you realize, that a lot of pixels has been "skipped". Plus looking closely on the edge of the last circle, there is clearly visible, that some dots are off the actual circumference - another result of the inaccuracy...

I could possibly use a formula (x - xorig)^2 + (y - yorig)^2 = r^2 to check every possible pixel in a rectangular area around the center, slightly bigger, than a diameter of the circle, if it does, or does't fall onto the cirle's circumference. But that would be very slow to repeat it all the time, as the circle grows.

Is there something, that could be done better? Could anyone help me to improve this? I don't insist on anything from my solution at all, and will accept any other solution, as long as it gives the desired results => let me read values of all (or the vast majority - 95%+) pixels on a circumference of a circle with given center and radius. The faster, the better...

Answer 1

1) Build a list of pixels of the smallest radius circumference. It is enough to keep the first octant (range 0..Pi/4 in the 1st quadrant of coordinate system) of circle, and get symmetric points with reflections. You can use, for example, Bresenham circle algorithm or just circle equation.

2) For the next iteration walk through all coordinates in the list (use right one, if there are two points with the same Y value) and check whether right neighbor (or two neighbors!) lies inside the next radius. For the last point check also top, right-top neighbor (at Pi/4 diagonal). Insert good neighbors (one or two) into the next coordinate list.

 Example for Y=5.
 R=8   X=5,6 //note that (5,5) point is not inside r=7 circle
 R=9   X=7
 R=10  X=8
 R=11  X=9
 R=12  X=10
 R=13  X=11,12 //!
 R=14  X=13

With this approach you will use all the pixels in maximal radius circle without gaps, and checking process for list generation is rather fast.

Edit: Code implements slightly another approach, it uses lower line pixel limit to built upper line.

It generates circles in given range, paints them to psychedelic colors. All math is in integers, no floats, no trigonometric functions! Pixels are used only for demonstration purposes.

procedure TForm1.Button16Click(Sender: TObject);

  procedure FillCircles(CX, CY, RMin, RMax: Integer);

    //control painting, slow due to Pixels using
    procedure PaintPixels(XX, YY, rad: Integer);
    var
      Color: TColor;
      r, g, B: Byte;
    begin
      g := (rad mod 16) * 16;
      r := (rad mod 7) * 42;
      B := (rad mod 11) * 25;
      Color := RGB(r, g, B);
     // Memo1.Lines.Add(Format('%d  %d  %d', [rad, XX, YY]));
      Canvas.Pixels[CX + XX, CY + YY] := Color;
      Canvas.Pixels[CX - YY, CY + XX] := Color;
      Canvas.Pixels[CX - XX, CY - YY] := Color;
      Canvas.Pixels[CX + YY, CY - XX] := Color;
      if XX <> YY then begin
        Canvas.Pixels[CX + YY, CY + XX] := Color;
        Canvas.Pixels[CX - XX, CY + YY] := Color;
        Canvas.Pixels[CX - YY, CY - XX] := Color;
        Canvas.Pixels[CX + XX, CY - YY] := Color;
      end;
    end;

  var
    Pts: array of array [0 .. 1] of Integer;
    iR, iY, SqD, SqrLast, SqrCurr, MX, LX, cnt: Integer;
  begin
    SetLength(Pts, RMax);

    for iR := RMin to RMax do begin
      SqrLast := Sqr(iR - 1) + 1;
      SqrCurr := Sqr(iR);
      LX := iR; // the most left X to check

      for iY := 0 to RMax do begin
        cnt := 0;
        Pts[iY, 1] := 0; // no second point at this Y-line
        for MX := LX to LX + 1 do begin
          SqD := MX * MX + iY * iY;
          if InRange(SqD, SqrLast, SqrCurr) then begin
            Pts[iY, cnt] := MX;
            Inc(cnt);
          end;
        end;

        PaintPixels(Pts[iY, 0], iY, iR);
        if cnt = 2 then
          PaintPixels(Pts[iY, 1], iY, iR);

        LX := Pts[iY, 0] - 1; // update left limit
        if LX < iY then // angle Pi/4 is reached
          Break;
      end;
    end;
    // here Pts contains all point coordinates for current iR radius
    //if list is not needed, remove Pts, just use PaintPixels-like output
  end;

begin
  FillCircles(100, 100, 10, 100);
  //enlarge your first quadrant to check for missed points
  StretchBlt(Canvas.Handle, 0, 200, 800, 800, Canvas.Handle, 100, 100, 100,
    100, SRCCOPY);
end;

Answer 2

If you want to make your code faster, don't call trigonometric functions inside the inner loop, increment sin(angle) and cos(angle) using

sin(n*step)=sin((n-1)*step)*cos(step)+sin(step)*cos((n-1)*step)
cos(n*step)=cos((n-1)*step)*cos(step)-sin(step)*sin((n-1)*step)

that is

  ...
  for currentRadius:= 80 to maxRadius do
  begin

    sinangle:= 0;
    cosangle:= 1;
    step:= 1 / currentRadius; // ?
    sinstep:= sin(step);
    cosstep:= cos(step);
    while {? } do
    begin

      xI:= Round(centerX + currentRadius * cosangle);
      yI:= Round(centerY + currentRadius * sinangle);

      newsin:= sinangle*cosstep + sinstep*cosangle;
      newcos:= cosangle*cosstep - sinstep*sinangle;

      sinangle:= newsin;
      cosangle:= newcos;

      ...    
    end;

  end;

Answer 3

First of all: you want all the points on a circle circumference. If you use any (good) algorithm, also some built-in circle function, you get indeed all the points, since the circumference is connected. What your picture shows, there are holes between neighbour circles, say r=100 and r=101. This is so for circumference drawing functions.

Now if you want that the pixels in your pixel set to cover all the pixels with incrementing radii, you can simply use following approach:

1. Build a set of filled circle pixels, say r = 101
2. Build a set of filled circle pixel with r = 100
3. Exclude set 2 from set 1

Filled circle algorithm is generally more efficient and simpler than connected circumference so you'll not loose much performance.

So you get a circumference which is slightly thicker than 1 px, but this set will surely cover the surface with growing radii without any holes. But it can also happen that the set built in such a way has overlapping pixels with previous set (r-1), so you'd know it better if you test it.

PS: Also it is not clear how any trigonometric functions appear in your code. I don't know any effective circle algorithm which use anything other than square root.

Answer 4

Why don't you simply use more digits for Pi and stop rounding to improve accuracy?

Further I suggest you use subpixel coordinates to get more accurate intensity values if you can afford the interpolation.

It's also very uncommon to use degrees in calculations. I highly recommend using radians. Not sure which functions you use here but Delphi's cos and sin seem to expect radians!