# Fast way to get positions of “boxed” array elements?

Given a 2D array (arbitrary numeric / evaluates to numeric elements, simplified for examples, and can be arbitrary dimensions so long as the operation makes sense, I'm working with 2K X 2K in my application):

test={{2, 2, 2, 2, 1, 1, 1, 2, 2, 1}, {2, 2, 2, 2, 2, 1, 1, 1, 2, 2}, {2,
2, 2, 1, 1, 1, 1, 1, 1, 2}, {1, 2, 2, 1, 2, 1, 1, 1, 2, 1}, {1, 2,
1, 1, 1, 2, 2, 1, 2, 2}, {1, 1, 2, 1, 2, 1, 2, 1, 2, 1}, {1, 2, 1,
1, 1, 2, 2, 1, 2, 1}, {2, 1, 2, 2, 1, 2, 1, 2, 2, 2}, {2, 2, 1, 1,
1, 1, 1, 1, 2, 1}, {2, 1, 1, 2, 1, 1, 2, 1, 2, 2}};


in MatrixForm:

I need to find the position(s), if any, of elements that are "boxed" by the same value of the element, with corners/edges treated the traditional way (i.e., corners have just 3 neighbors, edges have 5).

So, in the example array, the result should be

{{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {3, 7}}


I'm using

boxed[array_] :=
SparseArray[Unitize[Subtract[MaxFilter[array, 1], MinFilter[array, 1]]],
Automatic, 1]["NonzeroPositions"]


which performs fairly well, curious if there's a more elegant/better performing method.

Update: Here are some timings of answers so far on the netbook (for comparison, times for what I've tested on the workstation are ~1/20 of these). Logarithmic scale, else two become noise:

The current scheme I'm fiddling with is ~twice as fast as my current posted above, so there's twenty minutes saved on app. runs, but I'm sure there's some more clever ideas yet to be seen.

The nice compiled recent answer is not suitable, since elements can (and often are) out of machine-precision range (and it would be nice to extend this to completely arbitrary elements, e.g. symbols, etc., but not required.)

-

Ended up with this, which fits my needs of handling numerics not of machine precision. A bit over twice as fast as my OP, so pleased camper:

boxxed[array_] :=
masks = {{1 ;; -3, 1 ;; -3}, {3 ;;, 3 ;;}, {3 ;;, 1 ;; -3}, {1 ;; -3, 3 ;;},
{2 ;; -2, 1 ;; -3}, {2 ;; -2, 3 ;;}, {1 ;; -3, 2 ;; -2}, {3 ;;, 2 ;; -2}}},
SparseArray[Unitize@Total@Unitize@Map[Subtract[array, w[[Sequence @@ #]]] &, masks],
Automatic, 1]["NonzeroPositions"]]


Thanks also to all for the answers!

-

I like your approach. This is probably not going to be faster, but is certainly different:

Position[ListConvolve[BoxMatrix[1], ArrayPad[test, 1, f], {-1, 1}, {},
Times, Union@Flatten@{##} &] /. f -> Sequence[], {_Integer}]
(* {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {2, 2}, {3, 7}} *)


You can localize f if necessary.

-
+1 for convolution use, I used that in initial forays, and found (as you know I'm sure) that ListConvolve gets deadly slow when "custom" operators are used, as in 2+ orders of magnitude in this use. – ciao Mar 7 '14 at 4:26
@rasher one of the disappointing things about both ListConvolve and ListCorrelate is they are effectively unusable when you use custom functions. – Mike Honeychurch Mar 7 '14 at 7:53
@MikeHoneychurch: Yep, I've wondered why it takes such a gruesome hit, because even if the hit were 100% it could be hugely useful. – ciao Mar 7 '14 at 11:26
@rasher couldn't agree more. – Mike Honeychurch Mar 7 '14 at 21:28
@rasher I agree. For this particular case, if you leave out the corners/edges and only look at full 8-neighboured elements, then it's possible to do it without custom functions. – R. M. Mar 7 '14 at 21:31

I think this is one of the cases where a brute force compile-to-C might well win in the end, especially as you can reduce the number of copies with it. Here is one very straightforward try which seems to be quite fast already:

boxedC = Compile[{{array, _Integer, 2}},
Module[{
res, nx = Length[array], ny = Length[array[[1]]], im1, ip1, jm1, jp1
},
res = Table[0, {nx}, {ny}];
Do[
im1 = Max[i - 1, 1];
ip1 = Min[i + 1, nx];
Do[
jm1 = Max[j - 1, 1];
jp1 = Min[j + 1, ny];
If[And[
CompileGetElement[array, i, j] == CompileGetElement[array, im1, jm1],
CompileGetElement[array, i, j] == CompileGetElement[array, im1, j],
CompileGetElement[array, i, j] == CompileGetElement[array, im1, jp1],
CompileGetElement[array, i, j] == CompileGetElement[array, i, jm1],
CompileGetElement[array, i, j] == CompileGetElement[array, i, jp1],
CompileGetElement[array, i, j] == CompileGetElement[array, ip1, jm1],
CompileGetElement[array, i, j] == CompileGetElement[array, ip1, j],
CompileGetElement[array, i, j] == CompileGetElement[array, ip1, jp1]
],
res[[i, j]] = 1;
],
{j, 1, ny}
],
{i, 1, nx}
];
Position[res, 1]
],
CompilationTarget -> "C", RuntimeOptions -> "Speed"
];


One could probably earn some extra cycles with a smarter way to make those comparisons and handle the borders/corners in a more efficient way. It probably isn't worth the effort, though. And yes, it's not a very clever approach, especially as it copies the more subtle details from Leonid...

-
+1. I thought about taking this road too, but you beat me to it. Now it seems like we do have hear a full spectrum of solutions, for this problem. – Leonid Shifrin Mar 7 '14 at 15:17

Here's a ham-fisted approach...

Reap[
Do[
If[Equal @@
Flatten@test[[Max[1, i - 1] ;; Min[10, i + 1],
Max[1, j - 1] ;; Min[10, j + 1]]], Sow[{i, j}]];,
{i, 1, 10}, {j, 1, 10}]][[2, 1]]


Not too fast, not too pretty. Ham-fisted, indeed.

As rasher pointed out, method is hard-coded to input array size, it's easily modified to accept any size.

-
When I first ran this, I thought OMG... then realized you'd hard-coded array size. After adjustment, it's quickest answer other than OP so far. Thanks for posting. – ciao Mar 7 '14 at 11:25

I can't think of a way to improve the performance of your own code, but here is a more concise refactoring:

f2[array_] :=
With[
masks = Most @ Tuples[{;; -3, 3 ;;, 2 ;; -2}, 2]},
Unitize @ Total @ Unitize[ array ~Subtract~ w[[##]] & @@@ masks ]
] // SparseArray[#, Automatic, 1]["NonzeroPositions"] &


Or in a different style:

f3[array_] :=
ArrayPad[array, 1, "Fixed"] /. w_ :> (array ~Subtract~ w[[##]] &) @@@
Most @ Tuples[{;; -3, 3 ;;, 2 ;; -2}, 2] // Unitize // Total // Unitize //
SparseArray[#, Automatic, 1]["NonzeroPositions"] &

-
Ahh, nice! +1, I'd pondered for a moment how to condense the masks, since no perf. benefit, dropped it. You are a master of the terse... – ciao Mar 9 '14 at 9:38
@rasher Thanks for the +1. I just corrected the f3 code; I had a ( in the wrong place. Just goes to show that terse can be confusing -- explicit brackets do have value. :-) – Mr.Wizard Mar 9 '14 at 9:45

No idea how slow this is, but just for added variety:

SparseArray[
ImageData@
ImageFilter[If[Min[#] == Max[#], 1, 0] &, Image@test,
1]]["NonzeroPositions"]

-

This is not fast at all, but is worth mentioning as probably the easiest to adapt to finding other patterns:

Position[Partition[ArrayPad[test, 1, "Fixed"], {3, 3}, 1],
{{x_, x_, x_}, {x_, x_, x_}, {x_, x_, x_}}]

-

Yet another approach inspired by @simon-woods's excellent answer to this question is using the built-in LaplacianFilter without the need to pad the array as a pre-processing step.

This means that the nine-point stencil for the Laplace operator will be zero for 3x3 subarrays with all elements the same. I am pretty sure one of the answers includes this as an implementation (probably @rm-rf's), and it's still slower than the one you posted but it really is compact:

mybox[array_] := Position[N@LaplacianFilter[array, 1], 0.]


this works out of the box for symbolic matrices:

tst = RandomChoice[{x, y, z}, {50, 200}];


which is not the case with your boxxed as it stands.

I used N to force floating point arithmetic once symbols are encountered but it looks that this is done automatically if a numeric value is inserted i.e. try on

tst2 = RandomChoice[{1 + 10^(-320), 1 + 10^(-321)}, {10, 20}];


and you'll see what I mean. In short, I don't think you could go beyond machine precision with this approach.

---EDIT---

Stupidly, I didn't check if the nine-point stencil for the Laplacian gives false positives (it does). So using the same idea but different filters, this is essentially a slower version of Aky's:

 mybox[array_] := With[{prec = 30},
Position[
N[(MaxFilter[array, 1] - MinFilter[array, 1]), prec], N[0, prec]
]
]


but works just as well on symbolic and now arbitrary precision arrays.

-
Neat, but appears it gets wrong results on numerics. +1 for symbolic use though! – ciao Mar 10 '14 at 23:57
Works on numerics iff only two possible element values, so if you took my simplified example as test case, would sneak by. Still a cool way... – ciao Mar 11 '14 at 0:06
oops, good catch. I edited but now it really is the more general version of Aky's – gpap Mar 11 '14 at 10:41