Finding local maxima in multidimensional array (efficiently)

[Julian Francis]

Dear all,

How do I find the local maxima in a 4 Dimensional array.

e.g.

randTable = Table[Random[], {a, 1, 40}, {b, 1, 40}, {c, 1, 40}, {d, 1,
40}];

in this case there are 80 neighbours for each point (except at the
edges)

This is my best attempt, NB I take the highest 0.1% to reduce the
computational load:

In[11]:= threshold = Quantile[Flatten[randTable], .999]

Out[11]= 0.998989

In[13]:= top = Position[randTable, x_ /; x > threshold];

In[5]:= arrayBoundsCheck[array_, index_] :=
If[index == {}, True,
First[index] > 0 && First[index] <= Length[array] &&
arrayBoundsCheck[array[[First[index]]], Rest[index]]]

In[6]:= maxKernel =
DeleteCases[
Flatten[Table[{h, w, y, x}, {h, -1, 1}, {w, -1, 1}, {y, -1,
1}, {x, -1, 1}], 3], {0, 0, 0, 0}];

In[8]:= localMax[line_, array_] :=
Extract[array, line] >
Max[Extract[array,
Select[Map[Function[s, s + line], maxKernel],
arrayBoundsCheck[array, #] &]]]

In[15]:= Timing[results = Select[top, localMax[#, randTable] &];]

Out[15]= {80.59, Null}

Taking the top 0.1% as I have done is just done to try to bring the
compute time down, but isn't necessary to the problem.
I'm really just interested in finding a list of all the local maxima.
If someone has a much better solution, I would love to see it.

Thanks very much,
Julian.

[Ray Koopman]

Here are two routines that find local maxima.

'localMaxima' is a recoded version of your algorithm.

localMaxima[data_,p_] := Module[{dims,threshold,maxKernel,localMaxQ},
dims = Dimensions@data; threshold = Quantile[Flatten@data,p];
maxKernel = Transpose@Delete[Tuples[{-1,0,1},4],41];
localMaxQ[index_] := Extract[data, index] > Max@Extract[data,
Select[Transpose[index + maxKernel],
(And@@LessEqual@@@Transpose@{{1,1,1,1},#,dims}&)]];
Select[Position[data, x_ /; x > threshold], localMaxQ]]

'localMaxima2' takes a different approach. There is no thresholding,
and the window size is variable. The second parameter, 'h',
controls the window size: j +/- h on each axis.
A local maximum must be in the center of a full window.
('localMaxima' uses h = 1 implicitly, and does not
require maxima to be in the centers of full windows.)

localMaxima2[data_,h_] := Module[{dims = Dimensions@data,
r = Range[-h,h], s = (1+(1+2h)^4)/2}, Reap[
Do[If[#[[s]] > Max@Delete[#,s], Sow@{j1,j2,j3,j4}]& @
Flatten@data[[j1+r,j2+r,j3+r,j4+r]],
{j1,1+h,dims[[1]]-h},
{j2,1+h,dims[[2]]-h},
{j3,1+h,dims[[3]]-h},
{j4,1+h,dims[[4]]-h}]][[2,1]]]

Dimensions[randata = With[{n = 40}, Table[Random[],{n},{n},{n},{n}]]]
{40,40,40,40}

Timing@Length[result = localMaxima[randata,.999]]
{18.53 Second, 2448}

Timing@Length[result1 = localMaxima2[randata,1]]
{80.35 Second, 25788}

Timing@Length[result2 = localMaxima2[randata,2]]
{84.68 Second, 2629}

Timing@Length[result3 = localMaxima2[randata,3]]
{360.38 Second, 559}

[Ray Koopman]

'localMaxima3' defines the window around each test point by ANDing
a hypersphere of squared radius 'rr' with the hypercube used by
localMaxima2: only points that are in both the sphere and the cube
are in the window.

The minimum distance of the test point from each edge of the data
is still h, regardless of rr.

If rr >= 4 h^2 then localMaxima3 will give the same result as
localMaxima2 but will take more time to do it.

localMaxima3[data_,h_,rr_] := Module[{dims,g,sel,s},
dims = Dimensions@data; g = Range[-h,h];
sel = Flatten@Position[Tuples[g,4],x_/;x.x <= rr,{1}];
s = (1 + Length@sel)/2; Reap[

Do[If[#[[s]] > Max@Delete[#,s], Sow@{j1,j2,j3,j4}]& @

Flatten[data[[j1+g,j2+g,j3+g,j4+g]]][[sel]],

{j1,1+h,dims[[1]]-h},
{j2,1+h,dims[[2]]-h},
{j3,1+h,dims[[3]]-h},
{j4,1+h,dims[[4]]-h}]][[2,1]]]

Dimensions[randata = With[{n = 30},Table[Random[],{n},{n},{n},{n}]]]
{30,30,30,30}

Timing@Length[result = localMaxima[randata,.999]]

{5.45 Second, 772}

Timing@Length[res1 = localMaxima2[randata,1 ]]
Timing@Length[res104 = localMaxima3[randata,1,4]]
res1 == res104

{23.06 Second, 7614}
{25.66 Second, 7614}
True

Timing@Length[res2 = localMaxima2[randata,2 ]]
Timing@Length[res216 = localMaxima3[randata,2,16]]
Timing@Length[res204 = localMaxima3[randata,2, 4]]
res2 == res216

{22.69 Second, 701}
{31.02 Second, 701}
{22.98 Second, 5161}
True

Timing@Length[res3 = localMaxima2[randata,3 ]]
Timing@Length[res336 = localMaxima3[randata,3,36]]
Timing@Length[res325 = localMaxima3[randata,3,25]]
Timing@Length[res316 = localMaxima3[randata,3,16]]
Timing@Length[res309 = localMaxima3[randata,3, 9]]
res3 == res336

{ 89.30 Second, 142}
{144.59 Second, 142}
{141.58 Second, 155}
{ 80.03 Second, 254}
{ 65.54 Second, 776}
True

Timing@Length[res4 = localMaxima2[randata,4 ]]
Timing@Length[res464 = localMaxima3[randata,4,64]]
Timing@Length[res436 = localMaxima3[randata,4,36]]
Timing@Length[res425 = localMaxima3[randata,4,25]]
Timing@Length[res416 = localMaxima3[randata,4,16]]
res4 == res464

{143.83 Second, 40}
{265.62 Second, 40}
{242.99 Second, 47}
{174.37 Second, 86}
{105.37 Second,196}
True

[Julian Francis]

Ray,

Thank you very much for your solution, I had not thought of using Sow
and Reap in this way. It looks as if giving the algorithm domain
knowledge that all your local maxima will be above a certain value,
then thresholding first can save quite a bit of time. It does seem
quite a computationally demanding calculation.

I've just tested your localMaxima version; it gives the the same
answer on the test data set I used as my version. It seems
considerably faster than my original version (as well as more
compact); My timings are between 5-10 times faster than my original
coding.

Thanks for that,
Julian.

[Carl K. Woll]

On 11/28/2010 5:55 AM, Ray Koopman wrote:
> On Nov 27, 12:48 am, Ray Koopman<koop...@sfu.ca> wrote:
>> On Nov 26, 2:52 am, Julian Francis<julian.w.fran...@gmail.com> wrote:
>>> Dear all,
>>>
>>> How do I find the local maxima in a 4 Dimensional array.
>>>
>>> e.g.
>>>
>>> randTable = Table[Random[], {a, 1, 40}, {b, 1, 40}, {c, 1, 40}, {d, 1,
>>> 40}];
>>>
>>> in this case there are 80 neighbours for each point (except at the
>>> edges)
>>>
>>> This is my best attempt, NB I take the highest 0.1% to reduce the
>>> computational load:
>>>
>>> In[11]:= threshold = Quantile[Flatten[randTable], .999]
>>>
>>> Out[11]= 0.998989
>>>
>>> In[13]:= top = Position[randTable, x_ /; x> threshold];
>>>
>>> In[5]:= arrayBoundsCheck[array_, index_] :=
>>> If[index == {}, True,

>>> First[index]> 0&& First[index]<= Length[array]&&

If one only cares about rectangular blocks, then an alternative is to
use Partition, or even better, Developer`PartitionMap. For example, the
following finds the max value of each hypercube of length 2:

blockMax[data_] := Developer`PartitionMap[Max, data, ConstantArray[3,
Length[Dimensions[data]]], 1]

The following finds out the locations where the center of each hypercube
is at least as large as the max:

lm[data_] := Position[UnitStep[ArrayPad[data, -1] - blockMax[data]], 1] + 1

Comparing:

In[25]:= n=40;
data=RandomReal[1,{n,n,n,n}];

In[27]:= r1=lm[data];//Timing
r2=localMaxima2[data,1];//Timing
r1===r2

Out[27]= {3.385,Null}
Out[28]= {37.206,Null}
Out[29]= True

The above works because UnitStep[0.] is 1. If local max means that the
center value is strictly greater than all neighbors, then the above
approach can return incorrect results.

Carl Woll
Wolfram Research

[Julian Francis]

Dear Carl,

Thanks very much for your solution. I am getting somewhat different
timings than you have. You have lm working about ten times as fast as
localmaxima2, whereas on my machine I have lm about twice as fast as
localmaxima2. It would appear that you're machine is somewhat faster
than mine, but that doesn't explain the difference in the ratios of
the times for the two functions (I get 55 secs for localmaxima2 and 28
secs for lm).
I'm not sure what the reason for this is. My machine is a dual core
intel, and I am on Mathematica v7.

But thank you for the algorithm (apologies for late reply) but this is
an interesting approach, and I might adapt that for some other
problems I have. (My interests are machine vision).

Thanks,
Julian.