6
$\begingroup$

I am trying to find a faster way to do this than what I have put together or by looping through a bunch of points.

I want to select all pairs of integer points ((x1,y1),(x2,y2)) such that:

  1. Point (x1,y1) is at a specified distance r1 from center point (xc,yc).
  2. Point (x2,y2) is farther away from the center point than the point (x1,y1) is.
  3. Points (x1,y1) and (x2,y2) are separated by a specified distance r2.
  4. These points will be referring to 512X512 images, so each x and y value is limited to be between 1 and 512. Also because of this, I am having distances r1 and r2 only be integers, and I round any distances calculated between two points to the nearest integer.

Conceptually, for each point (x1, y1) that is a distance r1 from the center point, I am selecting all points that are a distance r2 from (x1, y1) and farther away from the center point. In the end, I want a list of pairs of points ((x1,y1),(x2,y2)). Here is my function:

points[r1_, r2_] :=
 Select[{x1, y1, x2, y2} /.
   FindInstance[
    1 <= x1 <= 512 &&
     1 <= y1 <= 512 &&
     1 <= x2 <= 512 &&
     1 <= y2 <= 512 &&
     Round[Sqrt[(x1 - xc)^2 + (y1 - yc)^2]] == r1 &&
     r1 <= Round[Sqrt[(x2 - xc)^2 + (y2 - yc)^2]] <= (r1 + r2),
    {x1, y1, x2, y2}, Integers, 512^2],
  Round[Sqrt[(#[[1]] - #[[3]])^2 + (#[[2]] - #[[4]])^2]] == r2 &]

I put condition 3 in the Select function rather than in the FindInstance function because the FindInstance function couldn't handle it. 512^2 maximum number of instances is a lot larger than I need, but it seems like having that number be too large doesn't hurt anything.

My code works, it is just very slow. The smallest r1 value I am using is 197, and r2 starts at 0 (r1 and r2 are limited to integers as well). My center point for now is at (229,256). I am sure there is an easier way to do this that is also faster. Any suggestions?


UPDATE:

Thanks to Michael E2 and Henrik Schumacher for helping with a faster solution. Although the below answers were not what I was looking for specifically, I was able to use them to obtain what I was looking for. Here is what I have put together:

pointFind[r_, c_] := 
  ArrayPad[(SparseArray@DiskMatrix[r]*
       SparseArray@(1 - DiskMatrix[r - 1, 2*r + 1])) // Transpose, 
    Transpose@{c - r - 1, 512 - c - r}]["NonzeroPositions"];

points[r1_, r2_, 
   c_] := ({#, 
       pointFind[r2, #]} /. {x_?NumericQ, y_?NumericQ} /; 
        Round[Norm[{x, y} - c]] < r1 -> Nothing) & /@ 
   pointFind[r1, c];

Here is an example of the result:

test = points[20, 10, {100, 100}];
Show[
 ListPlot[test[[;; , 1]],
  PlotRange -> {{65, 125}, {65, 125}},
  PlotStyle -> Black],
 ListPlot[test[[1, 2]],
  PlotStyle -> Red],
 Graphics[{Red, Point[test[[1, 1]]]}], AspectRatio -> 1]

I have plotted all points a rounded distance r1 (20) from the center, and then I picked one such point and show all of the points a distance r2 (10) from the first point that are farther from the center.

enter image description here

This works much much faster than what I had before. Thanks for everything!

$\endgroup$
  • $\begingroup$ What do you mean by "Select all pairs of points" - what will be returned by the function points? A Region or a pair of points randomly selected to fit the criteria? $\endgroup$ – Jason B. Dec 18 '17 at 21:04
  • $\begingroup$ Right now it returns a list of lists {x1,y1,x2,y2} that meet the criteria. $\endgroup$ – Aaron Stevens Dec 18 '17 at 21:07
  • $\begingroup$ The points are not being randomly selected. Given an r1 and r2, find ALL pairs ((x1,y1),(x2,y2)) that satisfy the conditions. $\endgroup$ – Aaron Stevens Dec 18 '17 at 21:32
  • $\begingroup$ do you really mean to find points at an exact integer distance r1? That changes the problem a good bit if you want exact vs allowing some tolerance. $\endgroup$ – george2079 Dec 19 '17 at 15:58
  • $\begingroup$ I guess technically I am allowing ranges of r-.5 and r+.5 for an integer r. That is why I say in my post that any distance I am measuring I round to the nearest integer. I just want to keep everything in terms of integers since I am dealing with pixels in an image. I am putting everything in terms of "to the nearest pixel". I am aware that not all distances will be exact integers, which is why I am rounding these distances to the nearest integer. The rounding is what supplies the tolerance. $\endgroup$ – Aaron Stevens Dec 19 '17 at 16:02
6
$\begingroup$

[Update 2: Used Nearest to compute final points pts. Another significant speed up.]

The OP's "distance" seems a bit odd and does not obey the triangle inequality. It leads to some odd points. But maybe that's what is being explored. For r1 = 220, r2 = 130 it takes 27 sec. 9 sec. 2.2 sec. Some of the black (x1,y1) points have no neighbor closer than Sqrt[2] which leads to occasional holes in the red (x2,y2) points.

A not very interesting brute-force method. [Update: I was working out a faster way, and during that time Henrik Schumacher also had the idea of using DiskMatrix. Multiplying the sparse arrays is slightly faster than subtracting the arrays. It's much faster than my original.]

[Update 3: In response to a comment asking for the pairs, I included two Association structures, p1$p2 to map {x1, y1} to the associated points {x2, y2}, and the inverse relation p2$p1.]

points2[{r1_, r2_}, c_] := 
  ArrayPad[(SparseArray@DiskMatrix[r2] *
       (SparseArray[1 - DiskMatrix[r1 - 1, 2 r2 + 1]])) // Transpose, 
    Transpose@{c - r2 - 1, 512 - c - r2}]["NonzeroPositions"];

Module[{r1 = 20, r2 = 10, c = {100, 100}, nf, p2$for$p1},
   pts1 = points2[{r1, r1}, c];
   pts2 = points2[{r1 + 1, r1 + r2}, c];
   nf = Nearest[pts1];
   p2$for$p1 = MapThread[     (* gets the points {x2, y2} for each point {x1, y1} *)
     Complement, {nf[pts2, {All, r2 + 1/2}], nf[pts2, {All, r2 - 1/2}]}];

   pts = Pick[pts2, Length /@ p2$for$p1 // Unitize, 1];  (* red points below *)
   p2$p1 = AssociationThread[pts2 -> p2$for$p1];         (* {x2, y2} -> {x1, y1} *)
   p1$p2 = Merge[Join]@ KeyValueMap[                     (* {x1, y1} -> {x2, y2} *)
      Function[{p2, p1}, AssociationThread[p1, \[FormalP]] /. \[FormalP] -> p2], 
      p2$p1];
   ]; // AbsoluteTiming
(*  {0.030094, Null} -- was {0.007771, Null} for just pts *)

Graphics[With[{r1 = 20, r2 = 10},
  {
   Circle[{100, 100}, r1 - 1/2], Circle[{100, 100}, r1 + 1/2],
   Point@pts1,
   Red, Circle[{100, 100}, r1 + r2 + 1/2],
   Point[pts]
   }
  ], GridLines -> {Range@512, Range@512}]

Mathematica graphics

Example use of the associations:

pts2[[138]]              (* pick a point {x2, y2} *)
(*  {75, 105}  *)

mypts2 = {75, 105};
mypts1 = p2$p1[mypts2]   (* gets the points {x2, y2} corresponding to {75, 105} *) 
Round[Norm[# - {100, 100}]] & /@ mypts  (* check distance *)
(*
{{80, 96}, {83, 111}}
{20, 20}
*)

p1$p2[First@mypts]       (* gets points {x2, y2} associated to {80, 96} *)
(*
{{70, 95}, {70, 96}, {70, 97}, {70, 98}, {70, 99}, {71, 91}, {71, 92}, {71, 100},
 {71, 101}, {72, 90}, {72, 102}, {73, 89}, {73, 103}, {74, 88}, {74, 104}, {75, 87},
 {75, 105}, {76, 87}, {76, 105}, {77, 86}, {77, 106}, {78, 86}, {78, 106}, {79, 86},
 {79, 106}, {80, 86}, {80, 106}, {81, 86}, {82, 86}, {83, 86}, {84, 87}}
*)

Norm[N[First@mypts - #]] & /@ p1$p2[First@mypts] // Round (* check distances *)
(*
{10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10}
*)
$\endgroup$
  • $\begingroup$ Hey Michael, something is odd here. When I run your code the image produced is quite different from the one I see here and I don't know how to interpret it... $\endgroup$ – Henrik Schumacher Dec 19 '17 at 6:07
  • $\begingroup$ @HenrikSchumacher Thanks, Henrik. I forgot to update the definition of points2[]. $\endgroup$ – Michael E2 Dec 19 '17 at 10:22
  • $\begingroup$ The speed looks very nice, and the output image looks like what I was imagining when visualizing all such pairs! I am going to try to unpack this right now to understand it more, but is there a way to modify this to keep the pairs together? I need the pairs of points ((x1,y1),(x2,y2)) kept together. Obviously, there will be multiple pairs corresponding to the same (x1,y1), but I need to keep them all together. I will go through it myself to see if I can do it, but if you realize it faster then I can go with that. $\endgroup$ – Aaron Stevens Dec 19 '17 at 15:09
  • $\begingroup$ @AaronStevens Please see the update. I wasn't sure in what form (data structure) you wanted the pairs. If you just want the straight pairs, I think this works: Table[ArrayPad[List /@ Complement[nf[p, {All, r2 + 1/2}], nf[p, {All, r2 - 1/2}]], {{0}, {1, 0}}, {{p}}], {p, pts2}] $\endgroup$ – Michael E2 Dec 20 '17 at 21:23
3
$\begingroup$

Maybe DiskMatrix is of interest for you. This is not a complete answer, but I try to show how you might apply DiskMatrix to your problem. For example, the set of all points {x1,y1} as specified above can be obtained with the following function:

digitalCirclePoints[xc_, yc_, a_, b_, r_] := 
 Module[{xlist, ylist, ipos, pat},
  {xlist, ylist} = Transpose[SparseArray[
      Subtract[
       DiskMatrix[r, {2 r + 1, 2 r + 1}],
       DiskMatrix[r - 1, {2 r + 1, 2 r + 1}]
       ]]["NonzeroPositions"]
    ];
  xlist = Clip[xlist + xc - (r + 1), {1, a}, {0, 0}];
  ylist = Clip[ylist + yc - (r + 1), {1, b}, {0, 0}];
  ipos = Flatten[Position[Positive[xlist ylist], True, 1]];
  pat = Transpose[{xlist[[ipos]], ylist[[ipos]]}]
  ]

Here are some usage examples:

SeedRandom[1234];
a = 64;
b = 48;
g = GraphicsGrid@Table[
   xc = RandomInteger[{1, a}];
   yc = RandomInteger[{1, b}];
   r1 = RandomInteger[{1, Mean[{a, b}]}];
   ArrayPlot[
    Transpose@
     SparseArray[digitalCirclePoints[xc, yc, a, b, r1] -> 1, {a, b}]],
   {3}, {3}]

enter image description here

Analogously, you can find the points {x2,y2} belonging to each {x1,y1} in a similar way.

$\endgroup$
0
$\begingroup$
centerpt = {5, 9};
r1 = 3;
r2 = 1;
fullpairs = {};
Table[
  firstpt = RandomPoint[Circle[centerpt, r1]];
  secondpt = RandomPoint[Circle[firstpt, r2]];
  If[EuclideanDistance[centerpt, secondpt] > 
     EuclideanDistance[centerpt, firstpt], 
   AppendTo[fullpairs, {firstpt, secondpt}]],
  {10}];
fullpairs

5000 points takes about seven seconds.

$\endgroup$
  • $\begingroup$ he wants integer pair i think. it is written in the title $\endgroup$ – Alucard Dec 18 '17 at 23:05
  • $\begingroup$ There is no guarantee that integer-coordinate points will obey the stated distance requirements as trivially demonstrated if $r2 = .5$. $\endgroup$ – David G. Stork Dec 18 '17 at 23:06
  • $\begingroup$ maybe he wants only those integer pairs that verify the relation at different distance r1 and r2 $\endgroup$ – Alucard Dec 18 '17 at 23:11
  • $\begingroup$ r2 must be a integer, at least it clearly must be in the OP's code. $\endgroup$ – Michael E2 Dec 18 '17 at 23:16
0
$\begingroup$

i am not sure this is completely right tho

pc = {229, 256};
r1 = 197;
r2 = 2;
v = {x, y};
lista = Select[{x, y} /. 
    Solve[EuclideanDistance[v, pc] == r1, v, Integers], 
   1 < #[[1]] < 512 && 1 < #[[2]] < 512 &];

Partition[
 Flatten[Table[
   Select[{x, y} /. 
     Solve[EuclideanDistance[v, lista[[k]]] == r2, v, Integers], 
    EuclideanDistance[pc, #] > 
      EuclideanDistance[pc, lista[[k]]] &], {k, 1, Length[lista]}]],
  2]
$\endgroup$
  • $\begingroup$ I think you forgot to count backwards... $\endgroup$ – Daniel Lichtblau Dec 18 '17 at 23:43
  • $\begingroup$ @DanielLichtblau what do you mean? I did not get that. $\endgroup$ – Alucard Dec 18 '17 at 23:50
  • 1
    $\begingroup$ Just based on the name $\endgroup$ – Daniel Lichtblau Dec 19 '17 at 20:20
0
$\begingroup$

here is a way to do this using region tools:

cp = {0, 0}
r1 = 10
r2 = 3
tol = .5
annulus[cp_, ri_, ro_] := 
 RegionDifference[Disk[cp, ro], Disk[cp, ri]]
p1 = {x, y} /.   FindInstance[
     Element[{x, y}, annulus[cp, r1 - tol, r1 + tol]], {x, y}, 
     Integers, 100];
result = {#, {x, y}} /. First@FindInstance[Element[{x, y},
        RegionIntersection[
         annulus[#, r2 - tol, r2 + tol],
         annulus[cp, Norm[#], Norm[#] + r2 + 2 tol]]], {x, y}, 
       Integers] & /@ p1;

Graphics[{Arrow /@ result, {Dashed, Circle[cp, r1]}}]

enter image description here

Edit..just noticed you want all such pairs. I dont think we can guarantee finding all this way but we can ask FindInstance to look for lots..

p1 = {x, y} /.  FindInstance[
     Element[{x, y}, annulus[cp, r1 - tol, r1 + tol]], {x, y}, 
     Integers, 
    Ceiling@Area@annulus[{0, 0}, r1 - tol, r1 + tol]];
Length@p1
result = Flatten[
  Function[{pp}, {pp, {x, y} /. #} & /@ 
     FindInstance[Element[{x, y},
         annulus[pp, r2 - tol, r2 + tol]
         ~RegionIntersection~
         annulus[cp, Norm[pp], Norm[pp] + r2 + 2 tol]], {x, y},
        Integers, 
           Ceiling@Area@annulus[{0, 0}, r2 - tol, r2 + tol]] ] /@ p1, 1]

this took 20 seconds to find 488 pairs.

enter image description here

for comparison a brute force version that gives the same result in .1 second.

  p1 = Select[Flatten[Table[{i, j} + cp,
       {i, -Ceiling[r1 + 1], Ceiling[r1 + 1]}, 
      {j, -Ceiling[r1 + 1], Ceiling[r1 + 1]}], 1] , 
    r1 - tol < Norm[# - cp] < r1 + tol &];
  result2 = 
   Flatten[Function[{pp}, {pp, #} & /@           
       Select[Flatten[Table[{i, j} + pp, 
          {i, -Ceiling[r2 + 1], Ceiling[r2 + 1]},
          {j, -Ceiling[r2 + 1], Ceiling[r2 + 1]}], 1] , (
          (r2 - tol < Norm[# - pp] < r2 + tol) && (Norm[# - cp] > 
             Norm[pp - cp])) &]] /@ p1, 1];

result == result2

True

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.