Re: Finding local maxima in multidimensional array (efficiently)
- To: mathgroup at smc.vnet.net
- Subject: [mg114224] Re: Finding local maxima in multidimensional array (efficiently)
- From: Ray Koopman <koopman at sfu.ca>
- Date: Sat, 27 Nov 2010 03:48:05 -0500 (EST)
- References: <ico3gv$olk$1@smc.vnet.net>
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}