# Faster way to pick pairs of integer coordinates from a "region"

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[(#[] - #[])^2 + (#[] - #[])^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_] :=
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. This works much much faster than what I had before. Thanks for everything!

• 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? Dec 18 '17 at 21:04
• Right now it returns a list of lists {x1,y1,x2,y2} that meet the criteria. Dec 18 '17 at 21:07
• The points are not being randomly selected. Given an r1 and r2, find ALL pairs ((x1,y1),(x2,y2)) that satisfy the conditions. Dec 18 '17 at 21:32
• 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. Dec 19 '17 at 15:58
• 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. Dec 19 '17 at 16:02

[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 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_] :=
(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}] Example use of the associations: pts2[] (* 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}
*)

• 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... Dec 19 '17 at 6:07
• @HenrikSchumacher Thanks, Henrik. I forgot to update the definition of points2[]. Dec 19 '17 at 10:22
• 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. Dec 19 '17 at 15:09
• @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}] Dec 20 '17 at 21:23

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;
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}] Analogously, you can find the points {x2,y2} belonging to each {x1,y1} in a similar way.

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.

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

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 < #[] < 512 && 1 < #[] < 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]

• I think you forgot to count backwards... Dec 18 '17 at 23:43
• @DanielLichtblau what do you mean? I did not get that. Dec 18 '17 at 23:50
• Just based on the name Dec 19 '17 at 20:20

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]}}] 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. 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