Re: Finding local maxima in multidimensional array (efficiently)
- To: mathgroup at smc.vnet.net
- Subject: [mg114246] Re: Finding local maxima in multidimensional array (efficiently)
- From: Ray Koopman <koopman at sfu.ca>
- Date: Sun, 28 Nov 2010 06:55:30 -0500 (EST)
- References: <ico3gv$olk$1@smc.vnet.net> <icqgl2$mfv$1@smc.vnet.net>
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