Re: Finding local maxima in multidimensional array (efficiently)
- To: mathgroup at smc.vnet.net
- Subject: [mg114327] Re: Finding local maxima in multidimensional array (efficiently)
- From: "Carl K. Woll" <carlw at wolfram.com>
- Date: Wed, 1 Dec 2010 02:11:20 -0500 (EST)
On 11/28/2010 5:55 AM, Ray Koopman wrote:
> On Nov 27, 12:48 am, Ray Koopman<koop... at sfu.ca> wrote:
>> On Nov 26, 2:52 am, Julian Francis<julian.w.fran... at 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]&&
>>> 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.
>> 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}
> '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
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